#### Study the asymmetry of the edges in the specialty referral matrix. 

We are studying a small graph, where each node is a primary specialty (e.g. PREVENTATIVE MEDICINE).  Then, an "edge" from one specialty to another specialty represents a referral from a doctor of the first type, to a doctor of the second type.  In a sample of 100k edges, there were 74 unique primary specialties. To accumulate edge strength, I used "two logs" as was done in the zip code graph.  You can [see my code here](code/smallSpecialtyReferralGraph.R).  You can download the data from the following url.
```{r}
rm(list = ls())
library(Matrix)
library(igraph)
load(url("http://pages.stat.wisc.edu/~karlrohe/netsci/data/specialty.RData"))
A = as.matrix(A) # it used to be a sparse matrix "dgCMatrix", even though it isn't that sparse. 
rs = rowSums(A); cs = colSums(A)
plot(rs,cs, xlab= "out degree", ylab = "in degree", col = "white")
text(rs, cs, labels=colnames(A), cex = .5)
ors = order(rs)
cbind(colnames(A)[ors], round(rs[ors],1))
```

Are edges reciprocal? 
```{r}
plot(A,t(A), main = "non-reciprocated edges at low degrees.")

bA = matrix(0, nrow = nrow(A), ncol = ncol(A))
bA[A>.5] = 1
nonrecA = bA*(1-t(bA))
nonrecEdges = which(nonrecA == 1, arr.ind = T)
# the first column makes some referrals to the second column, but it is much less likely in the other direction.
cbind(colnames(A)[nonrecEdges[,1]], colnames(A)[nonrecEdges[,2]])
```

What is the page rank? Imagine a Markov chain on the specialties.  You get randomly referred around this network.  Where are you likely to land?
```{r}
g= graph.adjacency(A,mode= "directed", weighted = T)
pr = page.rank(g,directed = T)
sort(pr$vec)
```

Before fitting the "fullly asymmetric" Stochastic Co-Blockmodel, let's try to fit the standard SBM allowing for $B$ to be asymmetric.  Use the standard spectral clustering algorithm applied to $A + A^T$.

```{r}
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:k1]
  if(length(k)==0){ 
    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)
Z = model.matrix(~as.factor(spec$clust)-1)
B = t(Z)%*%A%*%Z # this is counting the edges between groups.
round(B,2)  
table(spec$clust)  # should adjust for group size!
# TODO: Make a function that computes B for an SBM and DC-SBM given partition and graph.
image(t(Z)%*%As%*%Z, col = grey(seq(1,0,len=100)))
B - t(B)
```

Because of some stupid jpeg convention, R draws matrices not how we envision them. 
```{r}
imageMat = function(M){
  n = nrow(M)
  return(t(M[n:1,]))
}
image(imageMat(B), col = grey(seq(1,0,len=100)))

```


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

```{r}
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 = 3
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 this network was highly symmetric (i.e. reciprocal) what would that table look like?
```{r}
nm = colnames(A)
nm[uclust == 1]
nm[vclust == 1]
```
Note that $u$ corresponds to sending patterns and $v$ corresponds to receiving patterns.

What is a better 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."
```{r}
plot(s$d)
plot(s$d[-1])
k = 8 # remember to include the first one!
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))))
# running k-means on u and v together makes the cluster labels for u comparable to the labels for v.
km = kmeans(rbind(un,vn), centers = k,nstart = 100) 
uclust = km$clust[1:n]
vclust = km$clust[-(1:n)]
table(uclust,vclust)  # this is how **nodes** move between clusters.  
```
If this network was highly symmetric (i.e. reciprocal) what would that table look like?
```{r}
nm = colnames(A)
nm[uclust == 1]
nm[vclust == 1]
nm[uclust == 2]
nm[vclust == 2]
```
The "$B$" matrix gives the probabilities.  
```{r}
Y = model.matrix(~as.factor(uclust)-1)
Z = model.matrix(~as.factor(vclust)-1)  # dig into this code in class.

image(imageMat(t(Y)%*%A%*%Z), col = grey(seq(1,0,len=100)))
nm[uclust == 6]
```

