### Illustrating Stochastic Blockmodels (SBM's). 

In the first uses of the SBM, the block memberships $z$ were known a priori.  Later, $z$ was estimated in the  "*a posteriori* SBM".  In my reading of the literature, Newman and Girvan really pushed this idea further by saying that the "blocks" are modules.  You can think of the $B$ matrix as an adjacency matrix for a "modular" network; here, each block becomes a "supernode".  

It is reasonable to model $B$ using known covariates as $z$.  Then consider $\hat B$ as a weighted adjacency matrix.  In fact, we have done this twice in class!  In the first, we used billing zip codes.  In the second, we used "primary specialty" of physicians. Perhaps, this is intuitively similar to [vertex contractions](http://mathworld.wolfram.com/VertexContraction.html) in graph theory.  However, I suspect that the motivations are vastly different.  By "collapsing" on known blocks (e.g. zip code or specialty), we vastly simplify the network, allowing us to highlight salient features (e.g. geography or macro hospital referral patterns).

Sometimes, the modular network is sufficiently small that visualizations become tractable.  

```{r}
rm(list = ls())
library(Matrix)
library(igraph)
load(url("http://pages.stat.wisc.edu/~karlrohe/netsci/data/specialty.RData"))
g = graph.adjacency(A, mode = "directed", weighted = T)
plot(g, vertex.label = NA)
core = graph.coreness(g)
sum(E(g)$weight>3)
# E(g)$weight = as.numeric(E(g)$weight>3)
g = delete.edges(g, which(E(g)$weight <3))
plot(g, vertex.label = NA)
plot(induced.subgraph(g,vids = core>30),vertex.label=NA)
```

What a hairball!  Clearly, we need a more refined analysis (or better plotting software).  

Note: in this graph $\hat B$ is a type of weighted adjacency matrix.  Because $A$ is asymmetric, so is $\hat B$.  So, this modular network is directed. 

We have already modularized the network, but we can modularize it again!  This time, we can try to estimate a partition $z$ of specialty types.  

Recall that with a directed network, we can use the standard spectral clustering algorithm applied to $A + A^T$.



```{r}
# this function presumes that A is small and dense.  As such, it computes the full svd.
spectralClustering = function(A, k = c()){
  D = Diagonal(n = nrow(A), 1/sqrt(degree(g) + mean(degree(g))))
  L = D%*%A%*%D
  if(length(k) ==1) X = svd(L)$u[,1:k]
  if(length(k)==0){  # this will not work in markdown, but you can try.
    s = svd(L)
    plot(s$d[-1])
    print("type desired k (remember to add one)")
    k = scan()
    
    X = s$u[,1:k]
    }
  
  X = t(apply(X, 1, function(x) return(x/sqrt(sum(x^2)))))
  return(kmeans(X, k, nstart = 20))
  }

As = A + t(A)
spec = spectralClustering(As, 8)
table(spec$clust) # a bit unbalanced!
```

Note:  What is a good way to choose $k$?  Look at the scree plot!  From [wikipedia](https://en.wikipedia.org/wiki/Exploratory_factor_analysis#Cattell.27s_.281966.29_scree_plot), "Compute the eigenvalues for the correlation matrix and plot the values from largest to smallest. Examine the graph to determine the last substantial drop in the magnitude of eigenvalues. The number of plotted points before the last drop is the number of factors to include in the model. This method has been criticized because of its subjective nature (i.e., there is no clear objective definition of what constitutes a substantial drop).[17] As this procedure is subjective, Courtney (2013) does not recommend it."  Here is the scree plot that provided $k=8$.

```{r}
D = Diagonal(n = nrow(A), 1/sqrt(degree(g) + mean(degree(g))))
s= svd(D%*%A%*%D)
plot(s$d)
plot(s$d[-1])
```



So, now we found a partition into $k=8$ sets.  We can interpret the resulting partition with a $k = 8$ node graph!  To do this, we need to estimate $B$ for the given $z$.  For a given $z$, estimating $B$ is straightforward,
$$\hat B_{uv} = sum_{i,j: z(i) = u, z(j) = v}A_{ij}.$$
This can be done quickly with some linear algebra:

```{r}
fitSBM = function(A, z, y = c()){
  # takes label vector z (and y for co-blockmodel) as input!  Does not estimate partition.
  # computes the empirical probability of an edge from block u to block v 
  Zn = scale(model.matrix(~as.factor(z)-1),center=F, scale = T)/sqrt(length(z)-1)
  if(length(y)==0) return(t(Zn)%*%A%*%Zn)
  if(length(y)>0){
    Yn = scale(model.matrix(~as.factor(y)-1),center=F, scale = T)/sqrt(length(y)-1)
    return(t(Yn)%*%A%*%Zn)
    }
}
Bhat = as.matrix(fitSBM(A, spec$clust))  # convert to matrix from a Matrix.
round(Bhat,1)  
```
Huh?!  How come it has values greater than one??  Remember that the specialty graph $A$ is weighted.  So, "$B$"" contains the "expected weights".  As such, this can contain values greater than one!  If we were writing down likelihoods and fitting MLE's, we would have to be much more careful about all of this.  A great advantage of spectral techniques is that you can blindly apply them; of course, a downside is that you might forget and the code won't throw an error.

That table is difficult to read. Perhaps we could display as an image...

```{r}
image(Bhat, col = grey(seq(1,0,len=100)))
```

Because of some stupid jpeg convention, R draws matrices not how we envision them!  Be careful!!!

```{r}
imageMat = function(M){
  n = nrow(M)
  return(t(M[n:1,]))
}
image(imageMat(Bhat), col = grey(seq(1,0,len=100)))
image(imageMat((Bhat - t(Bhat))^2), col = grey(seq(1,0,len=100)))  # inspect the "reciprocity" 
```

Even with 8 node graphs, it is really difficult to simply visualize the adjacency matrix!

```{r}
gs = graph.adjacency(Bhat,mode = "directed", weighted=T)
loc = layout.fruchterman.reingold(gs)
plot(gs, layout= loc, edge.width=(E(gs)$weight), vertex.label = 1:8, edge.curved=T)
plot(gs, layout= loc, edge.width=sqrt(E(gs)$weight), vertex.label = table(spec$clust), edge.curved=T)
```
This is clearly a dense graph.  It is a bit difficult to interpret.  Best to focus on a node and move out from there.  The second block is highly connected! Guess which specialties are in this block!
```{r}
nm = colnames(A)
g = graph.adjacency(A, mode = "directed", weighted = T)
clustlist = list()
for(j in 1:8) clustlist[[j]]= c(
  round(mean(degree(g,v = spec$clust == j,mode = "in"))),
  round(mean(degree(g,v = spec$clust == j,mode = "out"))),  
  nm[spec$clust == j])
clustlist
```

Before reading the next steps, you might read the [lecture notes on the eigendecomposition and the svd](eigenSVDlecture.pdf).  
The algorithm "di-sim" explored here is described in [this paper](http://arxiv.org/abs/1204.2296). 

A key idea in the di-sim algorithm is that there are **two** partitions of the nodes. One for sending and one for recieving.  Each node is assigned to one sending-block and one receiving-block. this is encapsulated in the Stochastic Co-Blockmodel, where there are two partitions $Y$ and $Z$ and $E(A) = YBZ^T$. Before, the SBM provides a way of computing super nodes in a "collapsed" network. The ScBM summarizes a directed graph with a *bipartite* graph ([click here for more about bipartite](.  Each node is assigned to a "sending" super node and a "receiving" super node.  Let's illustrate this with the specialty graph.

```{r}
rs = rowSums(A); cs = colSums(A)
Drow = Diagonal(n = nrow(A), 1/sqrt(rs + mean(rs)))
Dcol = Diagonal(n = ncol(A), 1/sqrt(cs + mean(cs)))
L = Drow%*%A%*%Dcol
s = svd(L)
plot(s$d[-1])
```

```{r}
k = 8
u = s$u[,1:k]; v = s$v[,1:k]
n = nrow(A)
un = t(apply(u,1, function(x) return(x/sqrt(sum(x^2)+1/n))))
vn = t(apply(v,1, function(x) return(x/sqrt(sum(x^2)+1/n))))
km = kmeans(rbind(un,vn), centers = k,nstart = 100)
uclust = km$clust[1:n]
vclust = km$clust[-(1:n)]
table(uclust,vclust)
```
If "sending stochastic equivalence" implied "receiving stochastic equivalence", then this matrix would have a strong diagonal.  Let's visualize it with a bi-partite graph. 

```{r}
Bhat = fitSBM(A, y = uclust, z = vclust)
Ab = diag(rep(0, 2*k))
for(i in 1:k) for(j in 1:k) Ab[i,k+j] = Bhat[i,j]
Ab = Ab + t(Ab)
gb = graph.adjacency(Ab, weighted = T)
V(gb)$type = c(rep(T,8), rep(F,8))
V(gb)$name = c(paste("s",1:8, sep=""),paste("r",1:8, sep=""))
plot(gb, layout= layout.bipartite, edge.width=(E(gb)$weight)^.7, edge.arrow.size = 0)
```

What are the specialties in these blocks?

```{r}
nm = colnames(A)
sendlist = list()
reclist = list()
for(j in 1:8){
  sendlist[[j]] = nm[uclust == j]
  reclist[[j]] = nm[vclust == j]
  }
sendlist
reclist
```

