Skip to contents

Stochastic Block Models with covariates on nodes

S. Donnet

2026-11-10

Source: vignettes/SBM_nodes_covariates.qmd

1. Introduction

In a standard Stochastic Block Model (SBM), the nodes are assigned to blocks independently, with the same prior probabilities π=(π1,…,πQ)\pi = (\pi_1, \dots, \pi_Q) for all nodes. When some characteristics of the nodes are known (age, size, species traits…), it is natural to let these characteristics drive the probability of belonging to each block.

This vignette shows how to:

  1. simulate and fit simple and bipartite SBM in which the block memberships depend on covariates on nodes;
  2. check that taking the covariates into account can improve the clustering, especially when the network is noisy;
  3. use the Integrated Classification Likelihood (ICL) to decide whether a covariate has an effect on the clustering.

The inference relies on a variational EM algorithm implemented in the package blockmodels (Leger 2016).

We use the Adjusted Rand Index (ARI) from package aricode to compare an estimated clustering with the true one (ARI equals 11 for a perfect recovery, and is close to 00 for a random clustering).

ARI <- aricode::ARI

2. Simple SBM with covariates on nodes

2.1 The model

Consider a network of nn nodes, described by the adjacency matrix Y=(Yij)1≤i,j≤nY = (Y_{ij})_{1 \le i,j \le n}. Each node ii is described by a vector of pp covariates xi=(xi1,…,xip)⊤∈ℝpx_i = (x_{i1}, \dots, x_{ip})^\top \in \mathbb{R}^p. The first covariate is the intercept: xi1=1x_{i1} = 1 for all ii.

Each node ii belongs to one block Zi∈{1,…,Q}Z_i \in \{1, \dots, Q\}. The memberships are independent and follow a multinomial logistic regression on the covariates: ℙ(Zi=q∣xi)=πq(xi)=exp⁡(xi⊤βq)∑l=1Qexp⁡(xi⊤βl),q=1,…,Q, \mathbb{P}(Z_i = q \mid x_i) = \pi_q(x_i) = \frac{\exp(x_i^\top \beta_q)}{\sum_{l=1}^{Q} \exp(x_i^\top \beta_l)}, \qquad q = 1, \dots, Q, where βq∈ℝp\beta_q \in \mathbb{R}^p. For identifiability, the last block is the reference: βQ=0\beta_Q = 0. The parameter is gathered in the p×Qp \times Q matrix B=(β1,…,βQ)B = (\beta_1, \dots, \beta_Q), whose last column is null. The coefficient βkq\beta_{kq} is the effect of covariate kk on the log-odds of belonging to block qq rather than to block QQ: log⁡πq(xi)πQ(xi)=xi⊤βq=∑k=1pxikβkq. \log \frac{\pi_q(x_i)}{\pi_Q(x_i)} = x_i^\top \beta_q = \sum_{k=1}^p x_{ik} \beta_{kq}.

Given the blocks, the edges are independent and their distribution only depends on the blocks of the nodes they involve. In the Bernoulli case: Yij∣Zi=q,Zj=l∼ℬern(αql). Y_{ij} \mid Z_i = q, Z_j = l \sim \mathcal{B}ern(\alpha_{ql}). Poisson and Gaussian emission distributions are also available, as in the standard SBM.

The case with no covariates

When the only covariate is the intercept (p=1p = 1, xi=1x_i = 1), πq(xi)\pi_q(x_i) no longer depends on ii and we recover the standard SBM with πq=exp⁡(β1q)/∑lexp⁡(β1l)\pi_q = \exp(\beta_{1q}) / \sum_l \exp(\beta_{1l}).

2.2 Model selection: the ICL

The parameters θ=(B,α)\theta = (B, \alpha) are estimated with a variational EM algorithm, for a range of numbers of blocks QQ. The number of blocks, and more generally the model, is selected with the Integrated Classification Likelihood (ICL). For a model with QQ blocks and pp covariates on nodes, the ICL computed by the package is ICL(Q,p)=𝔼Z∼pVar[log⁡p(Y,Z∣X;θ̂)]−12[Q2log(n(n−1))+p(Q−1)logn], \mathrm{ICL}(Q, p) = \mathbb{E}_{Z\sim p^{Var}}[\log p(Y, Z \mid X; \widehat{\theta})] - \frac{1}{2}\left[ Q^2 \log\big(n(n-1)\big) + p\,(Q-1) \log n \right], where:

  • 𝔼Z∼pVar[log⁡p(Y,Z∣X;θ̂)]\mathbb{E}_{Z\sim p^{Var}}[\log p(Y, Z \mid X; \widehat{\theta})] is the integrated complete-data log-likelihood where the latent variables ZZ are integrated against to variational distribution pVarp^{Var}.
  • log⁡p(Y,Z∣X;θ)=log⁡p(Y∣Z;α)⏟network+log⁡p(Z∣X;B)⏟memberships; \log p(Y, Z \mid X; \theta) = \underbrace{\log p(Y \mid Z; \alpha)}_{\text{network}} + \underbrace{\log p(Z \mid X; B)}_{\text{memberships}};
  • Q2Q^2 is the number of connectivity parameters for a directed Bernoulli SBM (would be Q(Q+1)/2Q(Q+1)/2 for an undirected one), penalized by the logarithm of the number of dyads log⁡(n(n−1))\log (n(n-1)) (log⁡(n(n−1))/2\log (n(n-1))/2 for an undirected one);
  • p(Q−1)p(Q-1) is the number of free parameters in BB, penalized by the logarithm of the number of nodes nn.

Without covariates, p=1p = 1 and we recover the usual ICL of the SBM. Adding covariates increases the penalty, so the ICL favors the model with covariates only if they increase 𝔼Z∼pVar[log⁡p(Y,Z∣X;θ̂)]\mathbb{E}_{Z\sim p^{Var}}[\log p(Y, Z \mid X; \widehat{\theta})] enough, mainly through the term 𝔼Z∼pVar[log⁡p(Z∣X;θ̂)]\mathbb{E}_{Z\sim p^{Var}}[\log p(Z \mid X; \widehat{\theta})]. The number of blocks is chosen within each model, then the best model is the one with the highest ICL.

2.3 Simulation of a directed network

We simulate a directed binary network with n=120n = 120 nodes, Q=4Q = 4 blocks and p=3p = 3 covariates (the intercept and two Gaussian covariates).

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)
#>      intercept          x1         x2
#> [1,]         1 -0.96193342 -1.4954196
#> [2,]         1 -0.29252572  0.3678378
#> [3,]         1  0.25878822  0.5171440
#> [4,]         1 -1.15213189 -0.4843355
#> [5,]         1  0.19578283  0.6748556
#> [6,]         1  0.03012394 -0.7624486

The p×Qp \times Q matrix BB of the effects of the covariates is defined with the last block as reference (null last column).

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
#>           [,1]  [,2] [,3] [,4]
#> intercept  1.0 -2.01    2    0
#> x1         0.4  5.00   -2    0
#> x2         9.0  2.30    3    0

For this trial, the connectivity matrix α\alpha is drawn at random. The network is then simulated with sampleSimpleSBM, passing the covariates and their effects through the arguments nodesCovariates and nodesCovariatesParam.

Note

blockProp is not given: the block proportions are driven by the covariates.

alpha <- matrix(rbeta(Q * Q, 1, 8), Q, Q)
round(alpha, 2)
#>      [,1] [,2] [,3] [,4]
#> [1,] 0.16 0.04 0.23 0.05
#> [2,] 0.07 0.17 0.01 0.12
#> [3,] 0.01 0.06 0.09 0.00
#> [4,] 0.02 0.05 0.31 0.03

mySampler <- sampleSimpleSBM(
  nbNodes              = n,
  connectParam         = list(mean = alpha),
  model                = "bernoulli",
  directed             = TRUE,
  nodesCovariates      = X,
  nodesCovariatesParam = B
)
Z <- mySampler$memberships
Y <- mySampler$networkData
table(Z)
#> Z
#>  1  2  3  4 
#> 47 20 28 25
plot(mySampler)

To assess the methods in a more difficult setting, we also build a noisy version of the network, in which each entry of YY is flipped with probability ε=0.1\varepsilon = 0.1. The noisy network is still an SBM, with connectivity (1−ε)αql+ε(1−αql)(1 - \varepsilon) \alpha_{ql} + \varepsilon (1 - \alpha_{ql}): the blocks are less contrasted, hence harder to recover from the network only.

epsilon <- 0.1
Noise <- matrix(rbinom(n * n, 1, epsilon), n, n)
Ynoisy <- Y * (1 - Noise) + (1 - Y) * Noise

2.4 Estimation by VEM

The SBM with covariates on nodes is fitted with estimateSimpleSBM and the argument nodesCovariates.

Intercept

The matrix of covariates must contain the column of ones if an intercept is wanted: it is not added automatically.

estimOptions <- list(verbosity = 0, plot = FALSE)

fitCovar <- estimateSimpleSBM(
  Y, model = "bernoulli", directed = TRUE,
  nodesCovariates = X,
  estimOptions    = estimOptions
)
fitCovar
#> Fit of a Simple Stochastic Block Model -- bernoulli variant
#> =====================================================================
#> Dimension = ( 120 ) - ( 4 ) blocks.  No  edge covariate(s). 3  node covariate(s). 
#> =====================================================================
#> * Useful fields 
#>   $nbNodes, $modelName, $dimLabels, $nbBlocks, $nbCovariates, $nbDyads
#>   $blockProp, $connectParam, $covarParam, $covarList, $covarEffect 
#>   $nbNodesCovariates, $nodesCovariates, $nodesCovariatesParam 
#>   $expectation, $indMemberships, $memberships 
#> * R6 and S3 methods 
#>   $rNetwork, $rMemberships, $rEdges, plot(obj), print(obj), coef(obj) 
#> * Additional fields
#>   $probMemberships, $loglik, $ICL, $storedModels, 
#> * Additional methods 
#>   predict, fitted, $setModel, $reorder

For comparison, we fit the standard SBM, without covariates.

fitNoCovar <- estimateSimpleSBM(
  Y, model = "bernoulli", directed = TRUE,
  estimOptions = estimOptions
)

As usual, the number of blocks is selected with the ICL:

c(noCovar = fitNoCovar$nbBlocks, covar = fitCovar$nbBlocks)
#> noCovar   covar 
#>       4       4

and the estimated clusterings can be compared to the true one:

table(estimated = fitNoCovar$memberships, true = Z)
#>          true
#> estimated  1  2  3  4
#>         1 47  0  0  1
#>         2  0  0  0 24
#>         3  0 20  0  0
#>         4  0  0 28  0
table(estimated = fitCovar$memberships, true = Z)
#>          true
#> estimated  1  2  3  4
#>         1 47  0  0  0
#>         2  0  0  0 25
#>         3  0 20  0  0
#>         4  0  0 28  0

For this particular case, the clusters are nearly perfectly recovered, up to a label switching of the blocks (ARI = 0.98 without covariates and 1.00 with covariates).

2.5 Estimated effects of the covariates

The estimated matrix B̂\widehat{B} is given by the field nodesCovariatesParam, with the last block as reference.

round(fitCovar$nodesCovariatesParam, 2)
#>       [,1]  [,2]  [,3] [,4]
#> [1,] -0.51 -2.16 -3.12    0
#> [2,]  3.12  1.46  7.39    0
#> [3,]  4.61 -3.27  0.35    0

Label switching

The blocks being defined up to a permutation, the columns of B̂\widehat{B} are not in the same order as those of the true BB (used to simulate the data), and the reference block may differ. To compare them, we match each true block to its estimated counterpart, reorder the columns of B̂\widehat{B} accordingly and take the same reference block as in the simulation.

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)
#>           [,1]  [,2]  [,3] [,4]
#> intercept 1.65 -0.96  2.16    0
#> x1        1.66  5.93 -1.46    0
#> x2        7.88  3.62  3.27    0
B
#>           [,1]  [,2] [,3] [,4]
#> intercept  1.0 -2.01    2    0
#> x1         0.4  5.00   -2    0
#> x2         9.0  2.30    3    0

We can now compare the resulting probabilities ℙ̂(Zi=q)=π̂q(xi)\widehat{\mathbb{P}}(Z_i=q) = \widehat{\pi}_q(x_i) with their true values πq(xi)\pi_q(x_i).

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])))

2.6 Are the covariates selected?

On this network, the clusterings obtained with and without covariates are almost identical: the network alone is enough to recover the blocks. We can check whether the ICL favors the model with covariates. To do so, we also fit a model with two covariates drawn independently of everything else, hence unrelated to the blocks.

set.seed(123)
Xuseless <- cbind(1, matrix(rnorm(2 * n), ncol = 2))
fitUseless <- estimateSimpleSBM(
  Y, model = "bernoulli", directed = TRUE,
  nodesCovariates = Xuseless,
  estimOptions    = estimOptions
)
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)
covariates nbBlocks ARI ICL
none 4 0.98 -4106.70
x1 + x2 4 1.00 -4014.78
two useless covariates 4 0.98 -4117.72

The covariates barely change the clustering. The ICL selects the model with x1x_1 and x2x_2 (difference of ICL with the model without covariates: 91.9), whereas the useless covariates are not selected (difference of ICL: -11.0). To understand why, recall that the ICL (see above) is based on log⁡p(Y,Z∣X)\log p(Y, Z \mid X), which contains the term log⁡p(Z∣X)\log p(Z \mid X): without covariates it is ∑i∑qτ̂iqVElog⁡π̂q\sum_i \sum_q \hat{\tau}^{VE}_{iq} \log \widehat{\pi}_{q}, with covariates it is ∑i∑qτ̂iqVElog⁡π̂q(xi)\sum_i \sum_q \hat{\tau}^{VE}_{iq} \log \widehat{\pi}_{q}(x_i).

ElogpZ <- c(
  none      = sum(fitNoCovar$probMemberships%*%log(fitNoCovar$blockProp)),
  "x1 + x2" = sum(fitCovar$probMemberships*log(softmax(X %*% fitCovar$nodesCovariatesParam)))
)
ElogpZ
#>       none    x1 + x2 
#> -159.48117  -53.72192

The blocks being well predicted by the covariates, this term increases by 105.8, much more than the additional penalty 12(p−1)(Q−1)log⁡n=\frac{1}{2}(p-1)(Q-1)\log n = 14.4 due to the (p−1)(Q−1)(p-1)(Q-1) additional parameters.

What the ICL selects

The ICL selects the covariates that predict the block memberships, even when they do not improve the clustering, as here. Conversely, covariates unrelated to the blocks are not selected. Comparing the ICL of models with and without covariates therefore assesses whether the covariates explain the memberships, not whether they change the clustering.

2.7 Gain in clustering

We now fit both models on the noisy network and compare the ARI of the four estimated clusterings. The aim is to illustrate the fact that, if the covariates drive the blocks, they can help to recover the topology of the network.

fitNoCovarNoisy <- estimateSimpleSBM(
  Ynoisy, model = "bernoulli", directed = TRUE,
  estimOptions = estimOptions
)
fitCovarNoisy <- estimateSimpleSBM(
  Ynoisy, model = "bernoulli", directed = TRUE,
  nodesCovariates = X,
  estimOptions    = estimOptions
)
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)
network model nbBlocks ARI
original without covariates 4 0.98
original with covariates 4 1.00
noisy without covariates 2 0.20
noisy with covariates 4 0.94

On the original network, the structure is strong enough to be recovered from the network alone, and both models give similar clusterings. On the noisy network, the standard SBM selects only 2 blocks (ARI = 0.20): the network alone does not carry enough information. The covariates bring the missing information: the SBM with covariates selects 4 blocks (ARI = 0.94).

2.8 Testing the effect of covariates with the ICL

We now illustrate the fact that the ICL defined above can also be used to select the covariates that are involved in the block assignation.

We compare four models on the noisy network: no covariate (intercept only), x1x_1 only, x2x_2 only, and both. We also add a covariate x3x_3, drawn independently of everything else, which has no effect on the clustering.

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)
  }
})
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)
covariates nbBlocks ICL ARI deltaICL
none 2 -6608.50 0.20 0.00
x1 3 -6594.15 0.48 14.35
x2 4 -6584.17 0.79 24.33
x1 + x2 4 -6556.12 0.94 52.38
x1 + x2 + x3 4 -6563.28 0.94 45.22

deltaICL is the difference of ICL with the model without covariates: a positive value means that the covariates are useful to explain the clustering. Here, each of x1x_1 and x2x_2 has an effect on its own, and the best model (highest ICL) is the one with both of them, which is the true model. Adding the irrelevant covariate x3x_3 decreases the ICL: its small gain in fit does not compensate for the Q−1Q - 1 additional parameters.

Note

The ICL is not a statistical test in the strict sense (it does not provide a pp-value), but a model selection criterion that compares models of different complexities fitted on the same data.

3. Bipartite SBM with covariates on nodes

The same model and inference can be considered on bipartite models when covariates are provided, describing one or the two types of nodes.

3.1 The model

A bipartite network links n1n_1 nodes in rows to n2n_2 nodes in columns, through an incidence matrix Y=(Yij)1≤i≤n1,1≤j≤n2Y = (Y_{ij})_{1 \le i \le n_1, 1 \le j \le n_2}. Each row node ii (resp. column node jj) belongs to a block Zi1∈{1,…,Q1}Z^1_i \in \{1, \dots, Q_1\} (resp. Zj2∈{1,…,Q2}Z^2_j \in \{1, \dots, Q_2\}). Rows and columns can have their own covariates, xi∈ℝp1x_i \in \mathbb{R}^{p_1} and wj∈ℝp2w_j \in \mathbb{R}^{p_2}, acting through two multinomial logistic regressions: ℙ(Zi1=q∣xi)=exp⁡(xi⊤βq)∑l=1Q1exp⁡(xi⊤βl),ℙ(Zj2=r∣wj)=exp⁡(wj⊤γr)∑s=1Q2exp⁡(wj⊤γs), \mathbb{P}(Z^1_i = q \mid x_i) = \frac{\exp(x_i^\top \beta_q)}{\sum_{l=1}^{Q_1} \exp(x_i^\top \beta_l)}, \qquad \mathbb{P}(Z^2_j = r \mid w_j) = \frac{\exp(w_j^\top \gamma_r)}{\sum_{s=1}^{Q_2} \exp(w_j^\top \gamma_s)}, with βQ1=0\beta_{Q_1} = 0 and γQ2=0\gamma_{Q_2} = 0. All the memberships are independent and, given the blocks, the entries of YY are independent with Yij∣Zi1=q,Zj2=r∼ℬ(αqr). Y_{ij} \mid Z^1_i = q, Z^2_j = r \sim \mathcal{B}(\alpha_{qr}).

3.2 Simulation

We simulate a bipartite network with n1=120n_1 = 120 rows in Q1=4Q_1 = 4 blocks (with the same covariates design as before) and n2=100n_2 = 100 columns in Q2=3Q_2 = 3 blocks, described by an intercept and one Gaussian covariate. In the bipartite case, the covariates and their effects are given as lists of two matrices (rows, then columns).

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)
#> $row
#> 
#>  1  2  3  4 
#> 44 17 32 27 
#> 
#> $col
#> 
#>  1  2  3 
#> 22 39 39
plot(myBipSampler)

3.3 Missing data

On the bipartite network, instead of flipping the entries, we degrade the network by removing some of them: each entry of YY is set to NA independently with probability ε\varepsilon. We consider four levels of missing data: ε=10%\varepsilon = 10\%, 30%30\%, 60%60\% and 80%80\%.

Missing data in SBM

Missing entries are handled natively by the inference: the likelihood log⁡p(Y∣Z;α)\log p(Y \mid Z; \alpha) only involves the observed entries. Since here the entries are missing at random, independently of the blocks and of the network, the observed entries still follow the same SBM: the missing data do not bias the estimation, but they reduce the information available to recover the blocks.

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)))
#>     10% NA     30% NA     60% NA     80% NA 
#> 0.09658333 0.29783333 0.60058333 0.79883333

We plot the partially observed matrix where the non observed edges are plotted in white.

plotMyMatrix(YbipNA[[3]],clustering=Zbip)

3.4 Estimation and gain in clustering

The bipartite SBM with covariates on nodes is fitted with estimateBipartiteSBM, giving the list of row and column covariates in nodesCovariates. We fit the models with and without covariates on the original network and on each of the networks with missing data.

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)
  )
})
fitBip$original$covar
#> Fit of a Bipartite Stochastic Block Model -- bernoulli variant
#> =====================================================================
#> Dimension = ( 120 100 ) - ( 4 3 ) blocks.  No  edge covariate(s). 3 row node covariate(s). 2 col node covariate(s). 
#> =====================================================================
#> * Useful fields 
#>   $nbNodes, $modelName, $dimLabels, $nbBlocks, $nbCovariates, $nbDyads
#>   $blockProp, $connectParam, $covarParam, $covarList, $covarEffect 
#>   $nbNodesCovariates, $nodesCovariates, $nodesCovariatesParam 
#>   $expectation, $indMemberships, $memberships 
#> * R6 and S3 methods 
#>   $rNetwork, $rMemberships, $rEdges, plot(obj), print(obj), coef(obj) 
#> * Additional fields
#>   $probMemberships, $loglik, $ICL, $storedModels, 
#> * Additional methods 
#>   predict, fitted, $setModel, $reorder

The estimated effects of the covariates are given, for rows and columns, by nodesCovariatesParam:

lapply(fitBip$original$covar$nodesCovariatesParam, round, 2)
#> $row
#>       [,1]  [,2]  [,3] [,4]
#> [1,] -1.70 -2.91 -1.16    0
#> [2,]  3.41  7.04  3.51    0
#> [3,]  9.30 -1.39 -2.88    0
#> 
#> $col
#>       [,1]  [,2] [,3]
#> [1,] -1.43 -0.27    0
#> [2,]  2.00 -0.81    0
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)
network model nbRowBlocks nbColBlocks ARIrow ARIcol deltaICL
original without covariates 4 3 0.97 1.00 0.00
original with covariates 4 3 1.00 1.00 115.81
10% NA without covariates 4 3 0.97 1.00 0.00
10% NA with covariates 4 3 0.98 1.00 112.13
30% NA without covariates 4 3 0.93 1.00 0.00
30% NA with covariates 4 3 0.98 1.00 113.56
60% NA without covariates 3 3 0.69 0.85 0.00
60% NA with covariates 4 3 0.97 0.91 80.39
80% NA without covariates 1 1 0.00 0.00 0.00
80% NA with covariates 2 2 0.46 0.52 34.53

deltaICL is the difference of ICL with the model without covariates fitted on the same network.

Note

The ICL can only be compared between models fitted on the same network: with more missing entries, the likelihood involves fewer terms, so the ICL is not comparable across levels of missing data.

On the original network, the structure is strong enough to be recovered from the network alone: both models find the Q1=4Q_1 = 4 row blocks and the Q2=3Q_2 = 3 column blocks with almost perfect clusterings. The covariates do not improve the clustering, but the ICL still selects them, since they predict the memberships.

The SBM without covariates is robust to moderate levels of missing data: with 10% and 30% of missing entries, the clusterings are nearly unchanged (row ARI = 0.97 and 0.93). With 60% of missing entries, the network alone is no longer informative enough: the standard SBM selects 3 row blocks (ARI = 0.69), whereas the SBM with covariates finds 4 row blocks (ARI = 0.97). With 80% of missing entries, the network alone does not reveal any structure (1 ×\times 1 block), while the SBM with covariates finds 2 ×\times 2 blocks (ARI = 0.46 and 0.52 for rows and columns). The more entries are missing, the more the information brought by the covariates matters. For all the networks, the ICL selects the model with covariates.

3.5 Covariates on one type of nodes only

Covariates can be given for only one type of nodes, by setting an empty matrix for the other one. The column blocks then follow a standard multinomial distribution with proportions π2=(π12,…,πQ22)\pi^2 = (\pi^2_1, \dots, \pi^2_{Q_2}).

We illustrate it on the network with 60% of missing entries, for which the covariates improve the clustering.

fitBipRowOnly <- estimateBipartiteSBM(
  YbipNA[[lab[3]]], model = "bernoulli", dimLabels = c("row", "col"),
  nodesCovariates = list(X1, matrix(0, 0, 0)),
  estimOptions    = estimOptions
)
fitBipRowOnly$nbBlocks
#> row col 
#>   4   3

The ICL can then be used, as in the simple case, to decide which covariates are useful:

iclBip <- c(none        = fitBip[[lab[3]]]$noCovar$ICL,
            rowsOnly    = fitBipRowOnly$ICL,
            rowsAndCols = fitBip[[lab[3]]]$covar$ICL)
iclBip
#>        none    rowsOnly rowsAndCols 
#>   -2366.072   -2299.963   -2285.681

Here the highest ICL is obtained with covariates on both rows and columns.

References

Leger, Jean-Benoist. 2016. “Blockmodels: A r-Package for Estimating in Latent Block Model and Stochastic Block Model, with Various Probability Functions, with or Without Covariates.” arXiv Preprint arXiv:1602.07587.