knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
This vignette illustrates the use of the estimateMultiplexSBM
function and the methods accompanying the R6 classes multiplexSBMfit
.
The only package required for the analysis is sbm.
library(sbm)
Multiplex networks refer to a collection of networks involving the same sets of nodes, each network corresponding to a given type of interaction. Each network may also be referred as a layer. Such a network is studied in @multiplexSBM where the nodes are researchers and two networks are provided:
* a network representing the relations of advices between researchers, * a network representing the exchanges between the laboratories the researchers belong to.
In @multiplexKefi, the nodes are species and each network corresponds to a type of ecological relation, namely trophic links, negative non trophic links and positive non trophic links.
We give here a few details on Stochastic Block Models (SBM) for multiplex networks.
Assume that one studies $n$ nodes and observes $L$ types of links, resulting into $L$ matrices $(X^\ell_{ij}){i,j = 1,\dots,n}$, $\ell = 1, \dots, L$. Block models assume that the nodes are divided into $K$ clusters, this clustering being encoded into the latent variables $(Z_i){i=1,\dots,n}$ such that $Z_i = k$ if node $i$ belongs to cluster $k$.
The relations between nodes $i$ and $j$ are driven by the clusters they belong to :
$$ (X^1_{ij}, \dots,X^L_{ij}) | Z_i = q, Z_j = r \sim_{ind} \mathcal{F}(\cdot;\theta_{qr})$$ where $\mathcal{F}$ is a $L$ dimensional probability distribution adapted to the data.
Examples
Note that these two models assume that the $L$ interactions are dependent conditionally to the clustering. The integration against the latent variables $Z$ (clustering) implies that the interactions between two pairs of nodes are also dependent.
A collection of models are implemented in our package sbm
. These models can be classified into two groups : either the layers are dependent conditionally to $Z$ or not.
We consider conditional dependence in a few numbers of models :
for $L=2$ Binary networks (Bernoulli)
For any number $L$ of layers with Gaussian multivariate distributions but restricted to $\Sigma_{qr} = \Sigma$ (same covariance in any blocks).
In the Binary case, the number of parameters of dependent layers explodes with the number of layers $L$ and the results become difficult to understand. As a consequence, we also implemented inference methods for independent layers conditionally to the clusterings for any number of layers.
In the case of independent layers conditionally to the clusters, we assume that for any $\ell,i,j$, $$ X^\ell_{ij} | Z_i = q, Z_j = r \sim_{ind} \mathcal{F}\ell(\cdot;\theta^\ell{qr})$$ where any $\mathcal{F}_\ell$ can be Bernoulli, Poisson, Gaussian or Zero Inflated Gaussian.
All the previous models can be extended to bipartite networks (i.e. when the nodes in row are different from the nodes in column). In that case, two clusterings are introduced. Assume that the row (respectively column) nodes are divided into $K_1$ (resp. $K_2$) clusters, and denote $(Z^1_i){i=1,\dots,n}$ (resp. $(Z^2_j){i=j,\dots,p}$ ) the row (resp. col) clusterings.
The relations between nodes $i$ and $j$ are driven by the clusters they belong to :
$$ (X^1_{ij}, \dots,X^L_{ij}) | Z^1_i = q, Z^2_j = r \sim_{ind} \mathcal{F}(\cdot;\theta_{qr})$$
The maximization of the likelihood is obtained via a Variational version of the Expectation- Maximization algorithm.
The number of blocks is chosen via a penalized likelihood criterion (ICL).
See @multiplexSBM for details.
Note that, in the case of dependent networks, the inference is performed via the blockmodels
package, while the GREMLINS
package is used in the case of independent networks. As a consequence, the estimOptions
arguments are slightly different.
The function mySampleMultiplexSBM
supplies a method to simulate multiplex networks. The argument dependent
is set to FALSE
by default (assuming independent layers conditionally to the clusterings).
We simulate hereafter $2$ simple (i.e. non bipartite) layers. The corresponding matrices can be plotted using the function plotMyMultiplexMatrix
.
nbLayers <- 2
Examples
Nnodes <- 40 blockProp <- c(.4,.6) set.seed(1) connectParam <- list(list(mean=matrix(rbeta(4,.5,.5),2,2)),list(mean=matrix(rexp(4,.5),2,2))) model <- c("bernoulli","poisson") type <- "directed" mySampleMultiplexSBM_PB <- sampleMultiplexSBM( nbNodes = Nnodes, blockProp = blockProp, nbLayers = nbLayers, connectParam = connectParam, model=model, dimLabels = c('Individuals'), # generic name of the nodes type=type, seed = 1) listSBM_PB <- mySampleMultiplexSBM_PB$listSBM names(listSBM_PB) <- c("Bernoulli","Poisson") plotMyMultiplexMatrix(listSBM_PB)
The resulting object is a list of SBM objects (of the SBM class). To build such an object from an observed matrix use the function defineSBM
.
listSBM_PB
blockProp <- list(c(0.3, 0.3, 0.4), c(0.5, 0.5)) Q <- sapply(blockProp, function(p) length(p)) nbNodes <- c(80, 30) connectParam <- list() connectParam$mu <- vector("list", nbLayers) connectParam$mu[[1]] <- matrix(0.1, Q[1], Q[2]) + matrix(c(1, 1, 1, 0, 1, 0), Q[1], Q[2]) connectParam$mu[[2]] <- matrix(-2, Q[1], Q[2]) + matrix(c(1, 3, 2, 1, 2, 3), Q[1], Q[2]) connectParam$Sigma <- matrix(c(2, 1, 0.1, 4), nbLayers, nbLayers) model <- rep("gaussian", 2) mySampleMultiplexSBM_GG <- sampleMultiplexSBM( nbNodes = nbNodes, blockProp = blockProp, nbLayers = nbLayers, connectParam = connectParam, model = model, type = "bipartite", dependent = TRUE, dimLabels = c('row', 'col'), seed = 1) listSBM_GG <- mySampleMultiplexSBM_GG$listSBM plotMyMultiplexMatrix(listSBM_GG)
## MultiplexSBM Bernoulli with dependence Q <- 2 set.seed(94) P00 <- matrix(runif(Q * Q), Q, Q) P10 <- matrix(runif(Q * Q), Q, Q) P01 <- matrix(runif(Q * Q), Q, Q) P11 <- matrix(runif(Q * Q), Q, Q) SumP <- P00 + P10 + P01 + P11 P00 <- P00/SumP P01 <- P01/SumP P10 <- P10/SumP P11 <- P11/SumP connectParam <- list() connectParam$prob00 <- P00 connectParam$prob01 <- P01 connectParam$prob10 <- P10 connectParam$prob11 <- P11 model <- rep("bernoulli", 2) type <- "directed" nbLayers <- 2 Nnodes <- 40 blockProp <- c(0.6, 0.4) mySampleMultiplexSBM <- sampleMultiplexSBM( nbNodes = Nnodes, blockProp = blockProp, nbLayers = nbLayers, connectParam = connectParam, model = model, type = type, dependent = TRUE, seed = 1) listSBM_BB <- mySampleMultiplexSBM$listSBM plotMyMultiplexMatrix(listSBM_BB)
We are now able to perform inference on the multiplex network (search of the "best" number of clusters in terms of ICL). To do so, we have to chose the model, i.e. the distribution used in each matrix and the dependence or independence between matrices conditionally to the clusters.
load(file = 'resVignetteSimuMultiplex.rda')
res_PB <- estimateMultiplexSBM(listSBM_PB) res_PB$storedModels
res_PB$storedModels
One can now plot the reorganized matrices and the predicted values.
plot(res_PB)
plot(res_PB, type = 'expected')
One can also compare the estimated clusters to the simulated clusters. We recover the clusters perfectly.
All <- plotAlluvial(list(simulated = mySampleMultiplexSBM_PB$memberships$Individuals, estim = res_PB$memberships$Individuals)) All
- 2 Gaussian bipartite dependent Layers with 3 row clusters and 3 col clusters.
res_GG <- estimateMultiplexSBM(listSBM_GG, dependent = TRUE, estimOptions = list(plot = FALSE, verbosity = 0 )) res_GG$storedModels
res_GG$storedModels
plot(res_GG)
res_BB <- estimateMultiplexSBM(listSBM_BB,dependent = TRUE, estimOptions = list(plot = FALSE, verbosity = 0 )) res_BB$storedModels
res_BB$storedModels
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.