---
title: "Introduction to SBM and LBM. Practice"
author: "Sophie Donnet & Sarah Ouadah"
date: "Juillet 2020"
output: html_document
editor_options: 
  chunk_output_type: console
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE)
```




# Dataset : fungus-tree interactions  

VACHER, C., PIOU, D. and DESPREZ-LOUSTAU, M.-L. (2008). Architecture of an antagonistic tree/fungus network: The asymmetric influence of past evolutionary history. PLoS ONE 3 1740.

```{r import dataset,  echo = TRUE}
rm(list=ls())
load('inst/fungi_tree_data.Rdata')
ls()
```


- $154$ fungi species,  $51$ tree species
- "fungi_list": list of fungal species names
- "tree_list" : list of trees species names
- "tree_bin" : binary tree-tree interactions ($1$ if hosting at least one common fungal species, $0$ otherwise)
- "tree" : weighted tree-tree interactions (number of common fungal species two tree
species host) 
- "ListVar" : covariates discribing pairs of trees (genetic, taxonomic and geographic distances)
- "fungi_tree" : binary fungus-tree interactions

# Useful packages

- R-package blockmodels to infer Stochastic Block Model [SBM] and Latent Block model [LBM] on ecological networks
```{r install GSBM,  echo = TRUE,eval = FALSE, message=FALSE}
library(devtools)
install_github("GrossSBM/sbm")
```
(skip updates of dependecies if you have troubles)
```{r library  GSBM,  echo = TRUE,message=FALSE}
library(sbm)
```





GSBM is a common interface for all the pacakges developed in our unity (and friends) and dealing with the topological inference for complex networks. We use unified notations. 


- We will also need the ggplot2, igraph and alluvial packages for the graphical representations.

```{r load pack,  echo = TRUE,message=FALSE}
library(ggplot2)
library(igraph)
library(alluvial)
source('inst/function_for_blockmodels.R')
```

<!-- - For this class, we wrote additional ad-hoc functions. They are in the following file:  -->
<!-- ```{r source function,  echo = TRUE, message=FALSE} -->
<!-- source('function_for_blockmodels.R') -->
<!-- ``` -->


# Stochastic Block Model for Binary Network

## Dataset :  binary tree-tree interactions

We consider that two trees are in interaction if they host at least one common fungal species.

- Plot of the adjacency matrix:

```{r  tree tree bin plot, echo=TRUE} 

plotMatrix(Mat = tree_bin,rowFG = 'tree', colFG  = 'tree')
```


## Model and inference

\begin{align*}
 (Z_i) \text{ i.i.d.} \qquad & Z_i \sim \mathcal{M}(1, \pi) \\
 (Y_{ij}) \text{ indep.} \mid (Z_i) \qquad & (Y_{ij} \mid Z_i=k, Z_j = \ell) \sim \mathcal{B}(\alpha_{k\ell})
\end{align*}


simpleSBM refers to standard networks (w.r.t. bipartite)


- Then we apply the inference method
 
```{r estimate SBM, echo=TRUE, eval = TRUE}
sbm.tree_bin <- estimateSimpleSBM(
  tree_bin,
  model = "bernoulli", directed = FALSE
)
```


![](../Plot/ICL_tree_bin.png)

- The various points correspond to several runs of the estimation (resulting into a model selection criteria value) for a given number of clusters (denoted $Q$ here).
- The GSBM returns the better model (number of blocks selected with respect to the ICL criteria).

```{r select SBM, echo=TRUE, eval = TRUE}
sbm.tree_bin$nbBlocks
``` 


- We can now extract the parameters of interest, namely (`sbm.tree_bin$connectParam$mean`, $\alpha$, `sbm.tree_bin$blockProp` $\pi$) and the nodes (here tree species) clustering $Z$ (`sbm.tree_bin$memberships`) 
```{r extract param SBM, echo=TRUE, eval = TRUE}
sbm.tree_bin$blockProp
sbm.tree_bin$connectParam$mean
sbm.tree_bin$memberships
```

- One can plot the reorganized matrix according to the clustering.

```{r plot org treeb bin,  echo=TRUE, eval = TRUE}
plot(sbm.tree_bin,type = 'data',ordered = TRUE)
``` 

- Finally one can plot the macroview of the network.

```{r plot BM network treebin,  echo=TRUE, eval = TRUE}
G <- graph_from_adjacency_matrix(sbm.tree_bin$connectParam$mean, mode = c("undirected"), weighted = TRUE, diag = TRUE)
plot.igraph(G,vertex.size=sbm.tree_bin$blockProp*100,edge.width=abs(E(G)$weight)*3,vertex.color=1:sbm.tree_bin$nbBlocks, layout=layout_nicely)
``` 


- The composition of the clusters/blocks are given by:  
```{r list names blocks,  echo=TRUE, eval = TRUE}
lapply(1:sbm.tree_bin$nbBlocks,function(q){tree_list[sbm.tree_bin$memberships == q]})
```


#  SBM for Counting Network

## Dataset :  weighted tree-tree interactions

In the weighted tree-tree interactions network, each edge is valued with the number of common fungal species two tree species can host.

- Plot of the adjacency matrix:

```{r  tree tree plot, echo=TRUE} 
plotMatrix(Mat = tree,rowFG = 'tree', colFG  = 'tree')
```

## Model and inference
We consider the SBM Poisson model.

\begin{align*}
 (Z_i) \text{ i.i.d.} \qquad & Z_i \sim \mathcal{M}(1, \pi) \\
 (Y_{ij}) \text{ indep.} \mid (Z_i) \qquad & (Y_{ij} \mid Z_i=k, Z_j = \ell) \sim \mathcal{P}(\exp(\alpha_{kl})) = \mathcal{P}(\lambda_{kl})
\end{align*}

```{r SBM Poisson, echo=TRUE, eval = TRUE}
sbm.tree <- estimateSimpleSBM(
  tree,
  model = "poisson", directed = FALSE
)
```

![](../Plot/ICL_tree_poisson.png)

- We obtain the best numbers of clusters (with respect to the ICL criteria)
 


```{r select SBM poisson, echo=TRUE, eval = TRUE}
sbm.tree$nbBlocks
``` 

- We can now extract the parameters of interest, namely ($\lambda$, $\pi$) and the clustering of the nodes. 
```{r extract param SBM poisson, echo=TRUE, eval = TRUE}
sbm.tree$blockProp
sbm.tree$connectParam$mean
sbm.tree$memberships
```

- One can plot the reorganized matrix according to the clustering.

```{r plot org tree ,  echo=TRUE, eval = TRUE}
plot(sbm.tree, type = 'data',ordered = TRUE)
``` 

- One can plot the macroview of the network.


```{r plot BM network tree poisson, echo=TRUE, eval = FALSE}
G <- graph_from_adjacency_matrix(sbm.tree$connectParam$mean, mode = c("undirected"), weighted = TRUE, diag = TRUE)
plot.igraph(G,vertex.size=sbm.tree$blockProp*100,edge.width=abs(E(G)$weight)*3,vertex.color=1:sbm.tree$nbBlocks, layout=layout_nicely)
``` 
 

- The composition of the clusters/blocks are given by :  
```{r list names blocks Poisson,  echo=TRUE, eval = TRUE}
lapply(1:sbm.tree$nbBlocks,function(q){tree_list[sbm.tree$memberships == q]})
```

We are interested in comparing the two cluterings. To do so we use the alluvial flow plots.

```{r alluvial, echo=TRUE,eval=TRUE}
A <- as.data.frame(table(sbm.tree_bin$memberships,sbm.tree$memberships))
colnames(A)=c('SBM Bern',"SBM Poisson","Freq")
w   <- which(A$Freq!=0)
A <- A[w,]
alluvial(A[,c(1,2)],freq=A$Freq)
```


# SBM for Counting Network with covariates


## Dataset :  weighted tree-tree interactions with covariates

We have on each pair of trees 3 covariates, namely the genetic distance, the taxonomic distance and the geographic distance. 

Each covariate has to be introduced as a  matrix: $X^k_{ij}$ corresponds to the value of the $k$-th covariate describing the couple $(i,j)$. 

## Model and inference

We can also use the Blockmodels package to estimate the parameters described previously. 

\begin{align*}
 (Z_i) \text{ i.i.d.} \qquad & Z_i \sim \mathcal{M}(1, \pi) \\
 (Y_{ij}) \text{ indep.} \mid (Z_i) \qquad & (Y_{ij} \mid Z_i=k, Z_j = \ell) \sim \mathcal{P}(\exp(\alpha_{kl} + x_{ij}^\intercal \theta)) = \mathcal{P}(\lambda_{kl}\exp(x_{ij}^\intercal \theta))
\end{align*}



```{r covar SBM,echo=TRUE,eval= TRUE}
sbm.cov <- estimateSimpleSBM(netMat =tree,model = 'poisson',directed =FALSE, covariates  = ListVar)
```

- We select the best number of clusters (with respect to the ICL criteria)

```{r select SBM covar, echo=TRUE, eval = TRUE}
Q <- sbm.cov$nbBlocks
Q
``` 

- We can now extract the parameters of interest, namely ($\lambda$, $\pi$) and the clustering of the nodes. 
```{r extract param SBM poisson covar, echo=TRUE, eval = TRUE}
sbm.cov$connnectParam
sbm.cov$blockProp
sbm.cov$memberships
```


```{r plot org tree  cov,  echo=TRUE, eval = TRUE}
plot(sbm.cov,type='data', ordered = TRUE)
``` 

```{r plot BM network tree cov,  echo=TRUE, eval = FALSE}
G <- graph_from_adjacency_matrix(paramEstimSBMPoissonCov$lambda, mode = c("undirected"), weighted = TRUE, diag = TRUE)
plot.igraph(G,vertex.size = paramEstimSBMPoissonCov$pi*100,edge.width=sqrt(abs(E(G)$weight)),vertex.color = 1:Q, layout = layout_nicely)
``` 

![](../Plot/BMnet_poissoncov.png)

- The composition of the clusters/blocks are given by:  
```{r list names blocks poisson cov,  echo=TRUE, eval = TRUE}
lapply(1:sbm.cov$nbBlocks,function(q){tree_list[sbm.cov$memberships == q]})
```

 
We are interested in comparing the three cluterings. To do so we use the alluvial flow plots

```{r alluvial cov, echo=TRUE,eval=TRUE}
B <- as.data.frame(table(sbm.tree_bin$memberships,sbm.tree$memberships,sbm.cov$memberships))
colnames(B) = c('SBM Bern',"SBM Poisson","SBM Poisson + Cov","Freq")
w   <- which(B$Freq != 0)
B <- B[w,]
alluvial(B[,c(1,2,3)],freq = B$Freq)
```


# Latent Block model 

## Dataset : fungus-tree interactions
In the fungi-tree interactions network, an edge indicates if a tree species host a fungal species. 

- Plot of the incidence matrix

```{r  tree fungis plot, echo=TRUE} 
plotMatrix(Mat = fungi_tree,rowFG = 'fungi', colFG  = 'tree')
```

We can know try to get a clustering of fungus and tree species conjointly. 

## Model and inference 

\begin{align*}
 (Z^R_i) \text{ i.i.d.} \qquad & Z^R_i \sim \mathcal{M}(1, \pi^R) \\
  (Z^C_i) \text{ i.i.d.} \qquad & Z^C_i \sim \mathcal{M}(1, \pi^C) \\

 (Y_{ij}) \text{ indep.} \mid (Z^R_i, Z^C_j) \qquad & (Y_{ij} \mid Z^R_i=k, Z^C_j = \ell) \sim \mathcal{B}(\alpha_{k\ell})
\end{align*}



```{r LBM,echo=TRUE,eval=FALSE}
lbm <- estimateBipartiteSBM(as.matrix(fungi_tree),model = 'bernoulli')
```
```{r LBM save,echo=FALSE,eval=FALSE}
save(lbm,file='res_lbm.Rdata')
``` 

```{r LBM load,echo=FALSE,eval=TRUE}
load(file='../results/res_lbm.Rdata')
``` 


![](../Plot/ICL_fungi_tree.png)

- We select the best number of clusters (with respect to the ICL criteria)
```{r select lBM covar, echo=TRUE, eval = TRUE}
lbm$nbBlocks
```

- We can extract the parameters of interest, namely ($\alpha$, $\pi$) and the clustering of the nodes. 
```{r extract param LBM, echo=TRUE, eval = TRUE}
lbm$blockProp
lbm$connectParam
lbm$memberships$row
lbm$memberships$col
```


```{r plot org  tree fungis,  echo=TRUE, eval = TRUE}
plot(lbm, type ='data',ordered= TRUE)
``` 

```{r plot BM network tree fungi,  echo=TRUE, eval = TRUE}
G <- graph_from_incidence_matrix(lbm$connectParam$mean, weighted = TRUE)
plot(G,vertex.size=c(lbm$blockProp$row*100, lbm$blockProp$row*100), vertex.shape=c("circle", "square")[V(G)$type +1], edge.width=abs(E(G)$weight*2),vertex.color=1:sum(lbm$nbBlocks), layout=layout.bipartite)
``` 

- Exercise : give the name list of the various clusters. 



