TheSuperNetwork1.Rmd

---
title: "SuperNetwork1"
output: html_document
date: "2023-06-26"
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
  
## Working with the Super-Network 
  
In this doc we will explore working with the super-network, S, 
the mother of all networks. (More formally, the largest connected sub-graph in the union of all gene networks from KEGG.)

In this revision, we strip out all the things we are not using. (Older versions can now be found in
the Archive folder.)

Load whatever libraries we need:
```{r, load_libraries}
library(igraph)
library(ggplot2)
library("RColorBrewer") 
library("png")
library("ggraph")
library("networkD3")
library("animation")
library("maps")
library("geosphere")
library("RColorBrewer")
library(future.apply)
library(parallel)
```

Load the version of S constructed by Bryan:
```{r, loadS}
#load("~/Dropbox/Mac/Desktop/Simple_SIF_network_reactome.RData")
load("Simple_SIF_network_reactome.RData")
```

The network is in a variable called g.subg. It is an 'igraph'
We'll rename it to S:
```{r, rename}
S <- g.subg
rm(g.subg)
#head(S)
```

The summary of S gives us a few useful bits of info:
```{r, summary}
summary(S)
```

DN means that it is a 'D'irected network in which th e nodes have 'N'ames.
There are 11870 nodes and 347970 edges



Here's how to generate a random graph
(taken from https://r.igraph.org/articles/igraph.html).
```{r, random}
set.seed(33)
rg1 <- sample_grg(50, 0.2,coords=TRUE)
summary(rg1)
vertex_attr_names(rg1)
plot(rg1,frame.width=3)
plot(rg1,frame.width=6,label.cex=0.1,layout=layout_nicely)
```

That example generates a geometric random graph: n points are chosen randomly and uniformly inside the unit square and pairs of points closer to each other than a predefined distance d are connected by an edge. 


Manipulating the vertex attributes:
```{r, atrributes}
xcoord <- get.vertex.attribute(rg1,'x')
ycoord <- get.vertex.attribute(rg1,'y')
coords <- cbind(xcoord,ycoord)
plot(rg1,vertex.label.cex=0.5,vertex.layout=5*coords, vertex.size=5, vertex.color="blue",margin=0,asp=0.5)
V(rg1)$x   # to access the x-coords

```



Let's aim to construct a sub-graph of all vertices within 2 steps of a given vertex 
We try it on the super-network S.

First we write a function to pull out a sub-graph of all vertices within 'NeighborDist'
of vertex FocalVertex (the index of the focal vertex) on graph 'Network'
```{r, NbrFn}
NeighborSubgraph <- function(Network,FocalVertex,NeighborDist)
{
  DM <- distances(Network,v=V(Network),to=V(Network))
  sub <- as.numeric(DM[FocalVertex,] <= NeighborDist)
  subIDs <- which(sub==1)
  V(Network)$MyID <- 1:vcount(Network)
  ThisSub <- subgraph(Network,vids=subIDs)
  return (ThisSub)
}
```

And a function to pick low degree vertices, which might be good focal vertices:
```{r, lowdegreefn}
LowDegreeVertices <- function(Network,MaxDegree)
{
  IDs <- which(igraph::degree(Network)<=MaxDegree)
  #IDs <- which(LDVs==1)
  return (IDs)
}
```

Test it:
```{r , subgraphtest1}
#S2 <- make_ring(20)
set.seed(33)
S3 <- sample_grg(50, 0.2,coords=TRUE)
plot(S3)
x <- components(S3)
cat("Number of connected components: ",length(x$csize),"\n")
LowDegVerts <- LowDegreeVertices(S3,3)
(LowDegVerts)
SubNet <- NeighborSubgraph(S3,LowDegVerts[1],2)
plot(SubNet,vertex.label=V(SubNet)$MyID,vertex.label.cex=0.8,edge.arrow.size=0.3)
SubNet <- NeighborSubgraph(S3,LowDegVerts[1],3)
plot(SubNet,vertex.label=V(SubNet)$MyID,vertex.label.cex=0.8,edge.arrow.size=0.3)
SubNet <- NeighborSubgraph(S3,LowDegVerts[1],4)
plot(SubNet,vertex.label=V(SubNet)$MyID,vertex.label.cex=0.8,edge.arrow.size=0.3)

SubNet <- NeighborSubgraph(S3,LowDegVerts[9],4)
plot(SubNet,vertex.label=V(SubNet)$MyID,vertex.label.cex=0.8,edge.arrow.size=0.3)

```


Next, let's try defining a subgraph in a different way, which may create sub-graphs that are a bit
more like a piece of S.
Specifically, we start with one vertex and then iterate the following until the graph is big enough
1. Pick a vertex in the sub-graph uniformly at random
2. Find a neighbor of that vertex: add it if it is not already in the graph; else go back to 1.
We include all edges from the original graph that connect any two vertices in our sub-graph.
(This is handled automatically by the library's sub-graph function.)


```{r, IterativeBuildOfSubgraph}
BuildSubGraph <- function(g, subgsize) {
  size <- 0
  # label the vertices with our own IDs, that can be copied over
  V(g)$MyID <- 1:vcount(g)
  
  # pick a random vertex to start with (this may go wrong if the graph isn't fully-connected)
  subg <- sample(1:vcount(g), 1)
  size <- length(subg)
  itcount<-1
  #xcat("\nIteration: ",itcount,"  Subgraph size=",size,"  vertices: ",subg)
  while  ((size < subgsize)&&(itcount<1e6))
  {
    itcount <- itcount+1
    # pick a random vertex in the graph currently
    if (length(subg)==1){
      v<- subg
    }
    else
    {
      v <- sample(subg,1)
    }
    # get the neighbors of v
    v1 <- neighbors(g, V(g)[v], mode="all")
    # sample one of those
    if (length(v1)>0){
      if (length(v1)==1){
        v2 <- v1
      }else{
         v2 <- sample(v1, 1)
      }
      # add it to subg
      if (!(v2 %in% subg)){
        #cat("\nFocal vertex: ",v,"  nbrs: ",v1,"  sampled: ",v2)
        subg <- c(subg, v2)
      } 
      size <- length(subg)
    }
    #cat("\nIteration: ",itcount,"  Subgraph size=",size,"  vertices: ",subg)
    
    # debugging:
    compcount <- count_components(subgraph(g,vids=subg))
    if (compcount>1){
      cat("\n***Error. Disconnected subgraph")
      tempsub <- subgraph(g,vids=subg)
      #plot(tempsub,vertex.label=V(tempsub)$MyID,vertex.label.cex=0.8,edge.arrow.size=0.3)
     return(tempsub)
    }
  }
  Sub <- subgraph(g,vids=subg)
  if (size==subgsize){
    return(Sub)
  }else{
    cat("\nSub-graph routine failed: too many iterations")
    return(NA)
  }
}      
```

Test that function:
```{r, subgraphtest2}
set.seed(33)
GraphSize <- 100
ConnectionProb <- 0.2
g <- sample_grg(GraphSize, ConnectionProb,coords=TRUE)
while (count_components(g)>1){
  g <- sample_grg(GraphSize, ConnectionProb,coords=TRUE)
}
plot(g,vertex.label=V(g)$MyID)

SubGraphSize <- 10
for (i in 1:5){
  NewSub <- BuildSubGraph(g,SubGraphSize)
  plot(NewSub,vertex.label=V(NewSub)$MyID,vertex.label.cex=0.8,edge.arrow.size=0.3)
  cat("\n#components: ",count_components(NewSub))
  if (count_components(NewSub)>1){i<-5}
}
```

These all look good at first impression, and appear to be connected (as they must be, if we have done it correctly). But we should test the above more carefully to make sure it is really working.


Here's a smaller network to practice with. It comes from Bryan Queme. It's based on a single Alzheimer's pathway and is stored as an igraph called "alz"
```{r, SmallerNetwork}
load("Interconversion_of_nucleotide_di_and_triphosphates.RData")
summary(alz)
plot(alz,vertex.label.cex=0.2,edge.arrow.size=0.3,vertex.size=2)
```

It contains 551 vertices and 3789 edges. It's hard to see, so let's calculate some summary statistics:

```{r, summarystats}
get_diameter(alz, directed=F)

deg_alz <- igraph::degree(alz)
# two ways:
# simple way
hist(deg_alz,
     xlab = "degree",
     ylab = "Frequency",
     main = "Histogram of alz node degrees, without adjusting breaks and ylim",
     col = "skyblue")

# a bit more granularity
alz_table <- table(deg_alz)
relafreq <- alz_table/sum(alz_table)
barplot(relafreq, xlab = "degree",
        ylab = "Relative frequencies",
        main = "Degree distribution of alz",
        ylim = range(pretty(c(0,relafreq))),
        col = "orange")
```
***
&nbsp;

### Methods - Simple example of Local Moran's I and permutation-based testing

Now we begin to implement some of the statistical methodology.

We will practice with small random graphs, in which each possible edge 
occurs with probability 'EdgeProb'...
```{r, globals}
SizeOfOurNetwork <- 20
EdgeProb <- 0.15
```

We build another random network of 20 genes/nodes
```{r network, echo=FALSE}
RandomEdges <- function (x){
  x <- runif(1)<EdgeProb
  return (x)
}

RemoveSelfLoops <- function(x){
  if (x[1]==x[2])  x[3] <- 0
  return (x)
}

BuildRandomNetwork <- function(NetSize){
  MyNodes <- seq(1,NetSize)

  x <- seq(1, NetSize)
  y <- x
  
  MyEdges <- expand.grid(x = x, y = y)
  # make it undirected
  MyEdges2 <- MyEdges[MyEdges[,1]<MyEdges[,2],]

  EdgePresent <- rep(0,length(MyEdges2[,1]))
  EdgePresent <- apply(as.matrix(EdgePresent),MARGIN=1,FUN=RandomEdges)

  MyEdges3 <- cbind(MyEdges2,EdgePresent)
  
  MyEdges4 <- MyEdges3[MyEdges3[,3]==1,]

  ThisNetwork <- graph_from_data_frame(d=as.data.frame(MyEdges4),vertices=as.data.frame(MyNodes), directed=F) 
  
  return (ThisNetwork)
}

set.seed(593)

oldnet <- BuildRandomNetwork(SizeOfOurNetwork)
class(oldnet)

# plot it
l <- layout_with_fr(oldnet)
plot(oldnet, edge.arrow.size=.2, edge.curved=0,
     vertex.color="red", vertex.frame.color="white",
     vertex.label.cex=.7,vertex.size=2+igraph::degree(oldnet)) 
```

We can add attributes to the network, vertices, or edges as follows. Later, we will want the attribute for each node to correspond to gene expression values, for example.
```{r attributes}
net <- oldnet
V(net)$MyAttribute <- runif(length(V(net)),0,1)
#vertex_attr(net)
plot(net, edge.arrow.size=.5, vertex.label.color="black", vertex.label.dist=1.5,
vertex.size=4+10*V(net)$MyAttribute)
```


***

&nbsp;

## Taking advantage of what is there already...

It is simple to calculate distances between nodes:
```{r distancecalc}
(DM <- distances(net,v=V(net),to=V(net)))
```


Now we write a function to calculate Moran's (global) I. Reminder: Moran's-I is defined as \[ I=\frac{N}{W}\frac{\sum_i \sum_j w_{ij} (x_i-\bar{x})(x_j-\bar{x})}{\sum_i(x_i-\bar{x})^2} \]
This function, and what follows, needs to be carefully debugged to make sure it is calculating correctly.
```{r, MoransI, echo=TRUE}

MoransI <- function(DistanceMx, NodeAttributes){
  NoOfNodes <- length(DistanceMx[1,])
  WeightSum <- 0
  MoranSum <- 0
  DenomSum <- 0
  AttributeMean <- mean(NodeAttributes)
  #cat("\nAttribute mean= ",AttributeMean)
  for (i in 1:NoOfNodes){
    # the next line calculates the denominator
    DenomSum <- DenomSum + (NodeAttributes[i]-AttributeMean) * (NodeAttributes[i]-AttributeMean) 
    # and now we calculate the other bits
    for (j in 1:NoOfNodes){
      if ( i != j){
        ThisDist <- DistanceMx[i,j]
        if (ThisDist == 1) # neighbors only
        {
          WeightSum <- WeightSum + ThisDist # we use a weight of 1 if they are neighbors and 0 otherwise
          MoranSum <- MoranSum + ThisDist * (NodeAttributes[i]-AttributeMean) * (NodeAttributes[j]-AttributeMean)
        }
      }
    }
  }
  MoransI <- NoOfNodes * MoranSum / ( DenomSum * WeightSum)
}

cat("\nMorans-I: ",MoransI(DM,V(net)$MyAttribute),"    expectation= ",-1/(length(V(net)$MyAttribute)-1))

```


Test it out on 10 random graphs...
```{r tests, eval=FALSE,  echo=FALSE}
for (k in 1:10){
  net <- sample_gnp(SizeOfOurNetwork, EdgeProb,directed=FALSE)
  for (i in 1:length(V(net)))
    V(net)$MyAttribute[i] <- runif(1)
  vertex_attr(net)
  DM <- distances(net,v=V(net),to=V(net))
  plot(net, edge.arrow.size=.5, vertex.label.color="black", vertex.label.dist=1.5,
  vertex.size=4+10*V(net)$MyAttribute)
  
  cat("\nMorans-I: ",MoransI(DM,V(net)$MyAttribute),"    expectation= ",-1/(length(V(net)$MyAttribute)-1))

}

```

It doesn't seem to be matching the expectation very well, so let's keep an eye on that.

***

&nbsp;


Assessing the null distribution for Moran's-I via permutation tests.
We permute the node labels across the graph and then calc. Moran's-I to give us the null distribution:
```{r NullforMoransI, null, echo=FALSE}
PermutationTest <- function(ntwk,HowManyPermutations,WhichMeasure)
{
  Results <- rep(-9,HowManyPermutations)
  for (i in 1:HowManyPermutations){
    V(ntwk)$MyAttribute <- sample(V(ntwk)$MyAttribute,size=length(V(ntwk)$MyAttribute),replace=FALSE)
    
    if (WhichMeasure == 1)  # global Moran's-I
    {
      Results[i] <-  MoransI(DM,V(ntwk)$MyAttribute)
    }else{
      cat("\nUndefined measure for permuation test. Exit.")
      break;
    }
  }
  return (Results)
}

Results <-PermutationTest(net,200,1)
hist(Results,breaks=20,col="grey",main="Null dist. Obs(red)  Exp(green)")
abline(v=MoransI(DM,V(net)$MyAttribute),col="red")
abline(v=-1/(length(V(net)$MyAttribute)-1),col="seagreen")  # the expected value for Morans-I

```
This looks like a better match with the expectation, which is reassuring. 


How to find immediate neighbors of a vertex  (there is probably a built-in igraph function that does this automatically)
```{r, find_neighbors, echo=TRUE}
FindNeighbors <- function(ntwk,focnode)+{
  nbrs<-NULL
  for (i in 1:length(V(ntwk))){
    if ((distances(ntwk,v=(V(ntwk)==focnode),to=(V(ntwk)==i))==1) & (i!=focnode)){
      nbrs <- c(nbrs,i)
    }
  }
  return(nbrs)
}
plot(net)
(FindNeighbors(net,5))
(FindNeighbors(net,1))
```

Smoothing the labels to make them correlated. Here we use a simple proof of 
principle scheme in which we generate the labels independently and then smooth 
them by taking a weighted average of each label and the label of its neighboring 
vertices. Later on, we will try something more formal.
```{r, smoother, eval=TRUE, echo=FALSE}
Smoother <- function(ntwk,weight){
  NewLabels<-rep(0,length(V(ntwk)))
  for (i in 1:length(V(ntwk))){
    naybrs <- FindNeighbors(ntwk,i)
    NewL <- V(ntwk)$MyAttribute[i]
    for (j in naybrs){
      NewL <- NewL + weight * V(ntwk)$MyAttribute[j]
    }
    NewLabels[i]= NewL/(1+weight*length(naybrs)) 
  }
  return (NewLabels)
}

plot(net, edge.arrow.size=.5, vertex.label.color="black", vertex.label.dist=1.5,
vertex.size=4+10*V(net)$MyAttribute, main="Uniform attributes")

V(net)$MyAttribute <- Smoother(net,0.5)

plot(net, edge.arrow.size=.5, vertex.label.color="black", vertex.label.dist=1.5,
vertex.size=4+10*V(net)$MyAttribute,main="smoothed attributes")
```

### Test case: Global Moran's $I$

Let's compare Moran's-I for random labels and smoother labels...
```{r, compare, echo=FALSE}
NTimes <- 100
RandomMoransI <- rep(0,NTimes)
SmoothMoransI <- rep(0,NTimes)
for (i in 1:NTimes){
  # generate a random graph with random labels
  #net <- BuildRandomNetwork(SizeOfOurNetwork)
  net <- sample_gnp(SizeOfOurNetwork, EdgeProb,directed=FALSE)
  #for (j in 1:length(V(net)))
    V(net)$MyAttribute <- rnorm(length(V(net)),0,1)
  
  # calculate distances between nodes
  DM <- distances(net,v=V(net),to=V(net))

  # calculate Moran's-I for this graph
  RandomMoransI[i] <- MoransI(DM,V(net)$MyAttribute)
  
  # now smooth it and recalculate Moran's-I
  V(net)$MyAttribute <- Smoother(net,1)
  SmoothMoransI[i] <- MoransI(DM,V(net)$MyAttribute)
}

# compare via a violin plot
RandomMoransI <- cbind(rep("random",NTimes),RandomMoransI)
SmoothMoransI <- cbind(rep("smooth",NTimes),SmoothMoransI)
I <- rbind(RandomMoransI,SmoothMoransI)
dfsm <- data.frame( "Smooth" = I[,1], "MoransI" = as.numeric(I[,2]))

hw_p <- ggplot(dfsm, aes(x = Smooth, y = MoransI)) +
    geom_violin() + 
    geom_dotplot(binaxis='y', stackdir='center', dotsize=0.5, bin.width=60) +
    ggtitle("Moran's-I as a function of whether or not we smooth the vertex attributes") 
(hw_p)
```

Now that we've corrected the code for Moran's-I, these look quite different 

Let's try another way of producing correlated node labels....

***

&nbsp;

## Correlated node labels - spatially autocorrelated attribute values

And now for a more formal addition of spatial correlation, using correlated normals, where the degree of correlation depends upon the network structure (as originally proposed by George VY).
So, we generate data $y = (I_n - \rho W)^{-1} \times \epsilon$, where $\epsilon \sim MVN(0,I_n)$.
Since $\epsilon \sim MVN(0,I_n)$ the independent version has attrributes with $Normal(0,1)$ distribution.
```{r, SARmodel, echo=TRUE}
rho <- 0.8  # The degree of correlation 
SARsmoother <- function(ntwk,rho)
{
  # form edge matrix
  W <- as_adjacency_matrix(net, type = c("both"), names=FALSE, sparse=FALSE)
  n <- length(W[1,])
  #diag(W) <- 0
  W2 <- W/rowSums(W, na.rm = TRUE)

  # if nodes are disconnected, we will get NaNs, so set those to 0.
  W2[!is.finite(W2)] <- 0
  
  # generate the attributes
  y <- solve(diag(n) - rho * W2) %*% rnorm(n)
  
  # Check whether we generate spatial autocorrelation?
  library(ape)
  sc <- Moran.I(as.vector(y), W)
  cat("\np=",sc$p.value)

  return (y)
}

attribs <- SARsmoother(net,0.5)
net$MyAttributes <- SARsmoother(net,0.5)
```


Now write some tests for the above, comparing them to models in which the attributes 
are independent Normal(0,1).
```{r, SARtest, echo=FALSE}
NT <- 100
set.seed(49)

start_time <- Sys.time()

Idx <- 1:NT 
RandomMoransI <- rep(0,NT)
SARSmoothMoransI <- rep(0,NT)

aaa<-lapply(Idx, function(x){
  # generate a random graph with random labels
  net <- sample_gnp(SizeOfOurNetwork, EdgeProb,directed=FALSE)
  
  # calculate distances between nodes
  DM <- distances(net,v=V(net),to=V(net))

  # generate independent vertex attributes from a MVNormal(0,1)
  V(net)$MyAttribute <- rnorm(length(V(net)),0,1)

    # calculate Moran's-I for this graph
  RandomMoransI[x] <<- MoransI(DM,V(net)$MyAttribute)

  # now generate spatial correlated vertex labels and recalculate Moran's-I
  V(net)$MyAttribute <- SARsmoother(net,0.5)
  SARSmoothMoransI[x] <<- MoransI(DM,V(net)$MyAttribute)
})

# compare via a violin plot
RandomMoransI <- cbind(rep("random",NTimes),RandomMoransI)
SARSmoothMoransI <- cbind(rep("smooth",NTimes),SARSmoothMoransI)
I <- rbind(RandomMoransI,SARSmoothMoransI)
dfsm <- data.frame( "Smooth" = I[,1], "MoransI" = as.numeric(I[,2]))
#df$dataN <- as.factor(df$dataN)

hw_p <- ggplot(dfsm, aes(x = Smooth, y = MoransI)) +
    geom_violin() + 
    #geom_dotplot(binaxis='y', stackdir='center', dotsize=0.5, bin.width=60) +
    ggtitle("Moran's-I as a function of whether or not we smooth using SAR") 
(hw_p)

end_time <- Sys.time()
cat("\n time taken= ", end_time - start_time)

```

These don't seem to show much difference! So let's try something else. 
We pick a focal node at random and label it and all it's neighbors with a 5,
whereas everything else is labeled rnorm(0,1). That should show a difference.
First define the function to create the labels:

```{r, Correlated, echo=FALSE}
SpatiallyCorrelatedLabels <- function(ntwk,CentralNode){
  # label CentralNode and all its neighbors 1, and everything else Normal(0,1).
  V(ntwk)$MyAttribute <- rnorm(length(V(ntwk)),0,1)
  V(ntwk)$MyAttribute[CentralNode] <- 1
  for (i in 1:length(V(ntwk))){
    if (distances(ntwk,CentralNode,i) == 1){
      V(ntwk)$MyAttribute[i] <- 5
    }
    }
  return (ntwk)
}
```

Now test it out:
```{r, Test2, echo=FALSE}
set.seed(23)
# Pick a focal node at random
FocalNode <- sample(1:length(V(net)),1)
FocalNode <- 20
cat("\nFocalNode: ",FocalNode)
newnet <- SpatiallyCorrelatedLabels(net,FocalNode)
plot(newnet, vertex.color=ifelse(V(newnet)==FocalNode,"blue","red")) #, vertex.label.dist=1.5,vertex.size=4+10*V(newnet)$MyAttribute)
  
cat("\nMorans-I: ",MoransI(distances(newnet,v=V(newnet),to=V(newnet)),V(newnet)$MyAttribute),"    expectation= ",-1/(length(V(newnet)$MyAttribute)-1))
  
# Does it look significant?
Results2 <-PermutationTest(newnet,100,1)
hist(Results2,breaks=20,col="grey")
abline(v=MoransI(DM,V(newnet)$MyAttribute),col="red")
cat("\nObserved Mortan-s-I value is ",MoransI(DM,V(newnet)$MyAttribute))
```

Here we seem to see a reasonable amount of difference, provided the focal node is well-connected.

Let's look at local measures, which is what we think we will actually need in the end:

&nbsp;

***


## Local Moran's-I (LISA) - test cases.

We suppose that each node has some (binary or continuous) annotation $x_i$, 
and standardize those values by setting $z_i=x_i-\bar{x}$. 
The LISA measure of local clustering for each node, $i$, is then defined as 
$I_i = z_i \sum_{j \in J_i} w_{ij}z_j.$ 

Here, $J_i$ is the set of neighbors of node $i$ (although the definition 
can be generalized in an obvious way), $w_{ij}$ is a weight that is used to 
characterize the distance between nodes. For example, the weight might measure 
the number of edges on the shortest path between nodes $i$ and $j$.

We'll work with a bigger test network here

```{r, LISA, echo=FALSE}
set.seed(3975)
net <- sample_gnp(50, 0.15,directed=FALSE)
FocalNode <- 1
net <- SpatiallyCorrelatedLabels(net,1)
plot(net, vertex.color=ifelse(V(net)==FocalNode,"blue","red"),vertex.size=8+2*V(net)$MyAttribute)

MyLISA <- function(ntwk){
  V(ntwk)$LISAstat <- rep(-9,length(V(ntwk)))
  V(ntwk)$StandardizedAttribute <- rep(-9,length(V(ntwk)))
  AttributeMean <- mean(V(ntwk)$MyAttribute)
  # Standardize node labels
  V(ntwk)$StandardizedAttribute <- V(ntwk)$MyAttribute - AttributeMean
  for (i in 1:length(V(ntwk))){
    L <- 0
    for (j in 1:length(V(ntwk))){
      if (distances(ntwk,v=V(ntwk)[i],to=V(ntwk)[j]) == 1){
        # they are neighbors
        L <- L + V(ntwk)$StandardizedAttribute[j]
      }
    }
    L <- L * V(ntwk)$StandardizedAttribute[i]
    V(ntwk)$LISAstat[i] <- L
  }
  return (ntwk)
}

# calculate Local Moran's-I for each node and plot
LISAnet <- MyLISA(net)
plot(LISAnet, edge.arrow.size=.5, vertex.label.color="black", vertex.label.dist=1.5,
vertex.size=1+1*(0.01-min(V(LISAnet)$LISAstat)+V(LISAnet)$LISAstat)) 
cat("\nLISA stats for unpermuted data: ",V(LISAnet)$LISAstat)


LISAval <- V(LISAnet)$LISAstat
index <- seq(from=1, to=length(LISAval))
dataN <- rep(1,length(LISAval))
z1 <- cbind(dataN,index,LISAval)

# permute the attributes and repeat
PermNtwk <- net
V(PermNtwk)$MyAttribute <- sample(V(PermNtwk)$MyAttribute,size=length(V(PermNtwk)$MyAttribute),replace=FALSE)
LISAnetPerm <- MyLISA(PermNtwk)
plot(LISAnetPerm, edge.arrow.size=.5, vertex.label.color="black", vertex.label.dist=1.5,
vertex.size=1+(0.01-min(V(LISAnetPerm)$LISAstat)+V(LISAnetPerm)$LISAstat))
cat("\nLISA stats for permuted data: ",V(LISAnetPerm)$LISAstat)


LISAval<- V(LISAnetPerm)$LISAstat
index <- seq(from=1, to=length(LISAval))
dataN <- rep(2,length(LISAval))
z2 <- cbind(dataN,index,LISAval)

PermNtwk2 <- net
V(PermNtwk2)$MyAttribute <- sample(V(PermNtwk2)$MyAttribute,size=length(V(PermNtwk2)$MyAttribute),replace=FALSE)
LISAnetPerm <- MyLISA(PermNtwk2)
plot(LISAnetPerm, edge.arrow.size=.5, vertex.label.color="black", vertex.label.dist=1.5,
vertex.size=1+(0.01-min(V(LISAnetPerm)$LISAstat)+V(LISAnetPerm)$LISAstat))
cat("\nLISA stats for permuted data: ",V(LISAnetPerm)$LISAstat)
LISAval <- V(LISAnetPerm)$LISAstat
index <- seq(from=1, to=length(LISAval))
dataN <- rep(3,length(LISAval))
z3 <- cbind(dataN,index,LISAval)


## violin plot
zz <- rbind(z1,z2,z3)
df <- as.data.frame(zz)
df$dataN <- as.factor(df$dataN)

hw_p <- ggplot(df, aes(x = dataN, y = LISAval))
hw_p +
  geom_violin() + 
  geom_dotplot(binaxis='y', stackdir='center', dotsize=0.5) +
  ggtitle("Left=base; right=two permutations") +
  geom_hline(yintercept=V(LISAnet)$LISAstat[FocalNode], linetype=2, color="red", size=1)
 
  #scale_fill_brewer(palette="Dark2")
#violinplot(data=df$LISAval)
```
We see a clear difference there. the node that is the center of the cluster has a much higher LISA stat value. Let's investigate further...

# Semi-complete graphs

We will build a good network on which to explore correlated node labels. We will use a set of complete graphs joined to each other by a single vertex.
```{r,   SemiComplete, echo=FALSE}
BuildSemiCompleteNetwork <- function(SizeOfNetwork,CompleteSubnetworkSize){
  
  if (SizeOfNetwork%%CompleteSubnetworkSize!=0){
    cat("\nProblem with arguments to BuildSemiCompleteNetwork()")
    break
  }
  
  MyNodes <- seq(1,SizeOfNetwork)
  NumberOfCompleteSubgraphs <- SizeOfNetwork/CompleteSubnetworkSize

  x <- seq(1, SizeOfNetwork)
  y <- x
  
  MyEdges <- expand.grid(x = x, y = y)
  # make it undirected and with no self-loops
  MyEdges2 <- MyEdges[MyEdges[,1]<MyEdges[,2],]

  EdgePresent <- rep(0,length(MyEdges2[,1]))
  #EdgePresent <- apply(as.matrix(EdgePresent),MARGIN=1,FUN=RandomEdges)

  for (RowCounter in 1:length(MyEdges2[,1])){
    if ( floor((MyEdges2[RowCounter,1]-1)/CompleteSubnetworkSize)==floor((MyEdges2[RowCounter,2]-1)/CompleteSubnetworkSize) ){
        EdgePresent[RowCounter] <- 1
    }else{
      EdgePresent[RowCounter] <- 0
    }
  }
  
  # finally, add a connection between each complete subgraph and the next
  for (RowCounter in 1:length(MyEdges2[,1])){
    if  ( (MyEdges2[RowCounter,1]%%CompleteSubnetworkSize==1) && (MyEdges2[RowCounter,2]==MyEdges2[RowCounter,1]+CompleteSubnetworkSize)){
        EdgePresent[RowCounter] <- 1      
    }
    
  }

  MyEdges3 <- cbind(MyEdges2,EdgePresent)
  MyEdges4 <- MyEdges3[MyEdges3[,3]==1,]

  ThisNetwork <- graph_from_data_frame(d=as.data.frame(MyEdges4),vertices=as.data.frame(MyNodes), directed=F) 
  
  return (ThisNetwork)
  
}
```


Add correlated node labels to such a network
```{r, CorrelatedSemiComplete, echo=FALSE}
SpatiallyCorrelatedLabels1 <- function(ntwk,SizeOfSemiCompleteComponents){

  # label the first complete subgraph as N(8,1), and everything else N(1,1).
  V(ntwk)$MyAttribute <- 0
  for (i in 1:SizeOfSemiCompleteComponents){
    V(ntwk)$MyAttribute[i] <- rnorm(1,8,0.1)
  } 
  for (i in (SizeOfSemiCompleteComponents+1):length(V(ntwk))){
      V(ntwk)$MyAttribute[i] <- rnorm(1,1,1)
  } 
  return (ntwk)
}


  
  #cat("\nMorans-I: ",MoransI(distances(newnet,v=V(newnet),to=V(newnet)),V(newnet)$MyAttribute),"    expectation= ",-1/(length(V(newnet)$MyAttribute)-1))
  
  # Does it look significant?
  #Results2 <-PermutationTest(newnet,500,1)
#hist(Results2,breaks=20,col="grey")
#abline(v=MoransI(DM,V(newnet)$MyAttribute),col="red")


```

Test...
```{r, echo=FALSE}
SubNetworkSize <- 8
net<-BuildSemiCompleteNetwork(40,SubNetworkSize)
DM <- distances(net,v=V(net),to=V(net))  # distance mx

plot(net, edge.arrow.size=.5,  vertex.label.dist=1.5)
net<-SpatiallyCorrelatedLabels1(net,SubNetworkSize)
V(net)$MyAttribute
plot(net, edge.arrow.size=.5,  vertex.label.dist=1.5, vertex.size=5+2*V(net)$MyAttribute,main="attributes")
LISAnet <- MyLISA(net)
plot(LISAnet, edge.arrow.size=.5, vertex.label.color="black", vertex.label.dist=1.5,
vertex.size=0.25*(0.01-min(V(LISAnet)$LISAstat)+V(LISAnet)$LISAstat),main="LISA values")

# Moran-s I, for which we need a distance mx
cat("\nMoran's-I:",MoransI(DM,V(net)$MyAttribute))

# permute the attributes and repeat
PermNtwk <- net
V(PermNtwk)$MyAttribute <- sample(V(PermNtwk)$MyAttribute,size=length(V(PermNtwk)$MyAttribute),replace=FALSE)
plot(PermNtwk, edge.arrow.size=.5,  vertex.label.dist=1.5, vertex.size=5+2*V(PermNtwk)$MyAttribute,main="Permuted Attributes")
LISAnetPerm <- MyLISA(PermNtwk)
plot(LISAnetPerm, edge.arrow.size=.5, vertex.label.color="black", vertex.label.dist=1.5,
vertex.size=0.25*(0.01-min(V(LISAnetPerm)$LISAstat)+V(LISAnetPerm)$LISAstat),main="Permutated LISA values")
cat("\nPermuted Moran's-I:",MoransI(DM,V(PermNtwk)$MyAttribute))

```


Test for spatial structure by looking at the maximum value of the LISA stat across nodes for each graph
```{r, echo=FALSE}
set.seed(67)
net<-BuildSemiCompleteNetwork(40,SubNetworkSize)
net<-SpatiallyCorrelatedLabels1(net,SubNetworkSize)
DM <- distances(net,v=V(net),to=V(net))  # distance mx
LISAnet <- MyLISA(net)
plot(LISAnet, edge.arrow.size=.5, vertex.label.color="black", vertex.label.dist=1.5,
vertex.size=0.25*(0.01-min(V(LISAnet)$LISAstat)+V(LISAnet)$LISAstat),main="LISA stats on original network")
ObservedMaxLISA <- max(V(LISAnet)$LISAstat)
ObservedMoransI <- MoransI(DM,V(net)$MyAttribute)
NR <- 100
PermedMaxLISA <- rep(0,NR)
PermedMoransI <- rep(0,NR)
for (i in 1:NR){
  PermNtwk <- net
  V(PermNtwk)$MyAttribute <- sample(V(PermNtwk)$MyAttribute,size=length(V(PermNtwk)$MyAttribute),replace=FALSE)
  #plot(PermNtwk, edge.arrow.size=.5,  vertex.label.dist=1.5, vertex.size=5+2*V(PermNtwk)$MyAttribute,main="Permuted Attributes")
  LISAnetPerm <- MyLISA(PermNtwk)
  PermedMaxLISA[i] <- max(V(LISAnetPerm)$LISAstat)
  PermedMoransI[i] <- MoransI(DM,V(PermNtwk)$MyAttribute)
}
hist(PermedMaxLISA,main="Null distn of max(LISA stats); observed value shown in red",breaks=50,xlim=c(0,ObservedMaxLISA+10))
abline(v=ObservedMaxLISA,col="red",lty=2)
hist(PermedMoransI,main="Null distn of Moran's-I; observed value shown in red",breaks=50,xlim=c(0,max(ObservedMoransI,max(PermedMoransI))))
abline(v=ObservedMoransI,col="red",lty=2)
```

That seems to work nicely.

Here's a version in which we permute everything except the focal node when calculating the p-value for that node. (For now, let's deprecate this and just use permutations of all nodes. It will save time for testing and may well be the right thing to do anyway.)

First we need a permutation function that keeps node i fixed
```{r   OldFixedPerm, echo=FALSE, eval=TRUE}

FixedPerm <- function(Vect,FixedIndex){
  V <- sample(Vect[-FixedIndex],length(Vect)-1,replace=FALSE)
  W <- append(V,Vect[FixedIndex],after=FixedIndex-1)
}

CopyNet <- function(ntwk){
  return (ntwk)
}

SpatiallyCorrelatedLabels2 <- function(ntwk,SizeOfSemiCompleteComponents,Mean1,Mean2){

  # label the first complete subgraph as N(Mean1,1), and everything else N(Mean2,1).
  V(ntwk)$MyAttribute <- 0
  for (i in 1:SizeOfSemiCompleteComponents){
    V(ntwk)$MyAttribute[i] <- rnorm(1,Mean1,1)
  } 
  for (i in (SizeOfSemiCompleteComponents+1):length(V(ntwk))){
      V(ntwk)$MyAttribute[i] <- rnorm(1,Mean2,1)
  } 
  return (ntwk)
}

# Local moran's I (LISA) calculates a statistic for each node in the graph, indicating whether that node appears
# to be in a region of local clustering for the attribute (Pnode label) of interest.
# The following function calacluates the LISA stat for 'IthNode'.
# The statistic begins by subtracting the mean attrubute value from all attributes

LisaNetForNodeI <- function(ntwk,IthNode){
  V(ntwk)$StandardizedAttribute <- rep(-9,length(V(ntwk)))
  AttributeMean <- mean(V(ntwk)$MyAttribute)
  # Standardize node labels - subtract the mean attrubute value from all attributes
  V(ntwk)$StandardizedAttribute <- V(ntwk)$MyAttribute - AttributeMean
  for (i in IthNode:IthNode){  # This loop can be removed. It's a legacy of code I borrowed from a Morans-I function
    L <- 0
    for (j in 1:length(V(ntwk))){
      if (distances(ntwk,v=V(ntwk)[i],to=V(ntwk)[j]) == 1){
        # they are neighbors - for now we only consider directly adjacent nodes
        L <- L + V(ntwk)$StandardizedAttribute[j] # there would also be a weight term here if we wanted a more general definition of 'neighbor'
      }
    }
    L <- L * V(ntwk)$StandardizedAttribute[i]
  }
  return (L)
}

set.seed(45)
SizeOfOurNetwork <- 20
SubNetworkSize <- 5
net<-BuildSemiCompleteNetwork(SizeOfOurNetwork,SubNetworkSize)
DM <- distances(net,v=V(net),to=V(net))  # distance mx
net<-SpatiallyCorrelatedLabels2(net,SubNetworkSize,2,1)
V(net)$MyAttribute

plot(net, edge.arrow.size=.5,  vertex.label.dist=1.5, vertex.size=5+2*V(net)$MyAttribute,main="original attributes")

#LISAnet <- MyLISA(net)
#plot(LISAnet, edge.arrow.size=.5, vertex.label.color="black", vertex.label.dist=1.5,
#vertex.size=0.25*(0.01-min(V(LISAnet)$LISAstat)+V(LISAnet)$LISAstat),main="LISA values")

# permute the network and re-evaluate
i <- 1
V(net)$PermedAttribute <- FixedPerm(V(net)$MyAttribute,i)
V(net)$PermedAttribute
#plot(net, edge.arrow.size=.5,  vertex.label.dist=1.5, vertex.size=5+2*V(net)$PermedAttribute,main=paste("permuted attributes with node",i," fixed"))
#LISAnet2 <- MyLISA(net)

net2 <- CopyNet(net)
plot(net2, edge.arrow.size=.5,  vertex.label.dist=1.5, vertex.size=5+2*V(net2)$MyAttribute,main="attributes of permuted network")

net2 <- CopyNet(net)
NPerms <- 20
PValue <- rep(0,SizeOfOurNetwork)
LStat <- rep(0,NPerms)
for (i in 1:SizeOfOurNetwork)
{
  cat("\nPerm",i)
  # Evaluate Lisa stat for node i
  for (j in 1:NPerms){
    V(net2)$MyAttribute <- FixedPerm(V(net2)$MyAttribute,i)
    LStat[j] <- LisaNetForNodeI(net2,i)
  }
  #boxplot(LStat)
  abline(h=V(LISAnet)$LISAstat[i],col="red")
  PValue[i] <- sum(LStat>V(LISAnet)$LISAstat[i])/NPerms
}
(PValue)
#boxplot(PValue[1:25])
#boxplot(PValue[26:200])
OverExp <- c(rep(0,25),rep(1,175))
PlotData <- as.data.frame(cbind(PValue,OverExp))

ggplot(PlotData, aes(x = factor(OverExp), y = PValue)) +
  #geom_boxplot()  +
  geom_violin() +
   geom_point() +
   xlab("Group") +
  ylab("Permutation test P-value") +
   theme_grey(base_size = 24, base_family = "Times")
  # + theme(axis.title.x = element_text(face = "italic", colour = "black"))
  #+ theme(axis.text.x = element_text(family = "Times", size = rel(5)))
  #geom_violin()
  #  + scale_y_log10()
```

Here's a new version of the permutation test, written by Mingzhi: 

```{r   NewFixedPerm, echo=FALSE, eval=TRUE}

FixedPerm <- function(Vect,FixedIndex){
  V <- sample(Vect[-FixedIndex],length(Vect)-1,replace=FALSE)
  W <- append(V,Vect[FixedIndex],after=FixedIndex-1)
}

CopyNet <- function(ntwk){
  return (ntwk)
}

SpatiallyCorrelatedLabels2 <- function(ntwk,SizeOfSemiCompleteComponents,Mean1,Mean2){

  # label the first complete subgraph as N(Mean1,1), and everything else N(Mean2,1).
  V(ntwk)$MyAttribute <- 0
  for (i in 1:SizeOfSemiCompleteComponents){
    V(ntwk)$MyAttribute[i] <- rnorm(1,Mean1,1)
  } 
  for (i in (SizeOfSemiCompleteComponents+1):length(V(ntwk))){
      V(ntwk)$MyAttribute[i] <- rnorm(1,Mean2,1)
  } 
  return (ntwk)
}

# Local moran's I (LISA) calculates a statistic for each node in the graph, indicating whether that node appears
# to be in a region of local clustering for the attribute (Pnode label) of interest.
# The following function calacluates the LISA stat for 'IthNode'.
# The statistic begins by subtracting the mean attrubute value from all attributes
NeighborsLists <- replicate(length(V(net2)), c(), simplify = FALSE)
for(i in 1:length(E(net2))){
  a=as.integer(get.edgelist(net2)[i,1])
  b=as.integer(get.edgelist(net2)[i,2])
  NeighborsLists[[a]]<-c(NeighborsLists[[a]],b)
  NeighborsLists[[b]]<-c(NeighborsLists[[b]],a)
}

LisaNetForNodeI <- function(ntwk,IthNode){
  V(ntwk)$StandardizedAttribute <- rep(-9,length(V(ntwk)))
  AttributeMean <- mean(V(ntwk)$MyAttribute)
  # Standardize node labels - subtract the mean attrubute value from all attributes
  V(ntwk)$StandardizedAttribute <- V(ntwk)$MyAttribute - AttributeMean
  L <- 0
  for (j in NeighborsLists[[IthNode]]){
      # they are neighbors - for now we only consider directly adjacent nodes
    L <- L + V(ntwk)$StandardizedAttribute[j] # there would also be a weight term here if we wanted a more general definition of 'neighbor'
  }
  L <- L * V(ntwk)$StandardizedAttribute[IthNode]
  
  return (L)
}

set.seed(45)
SizeOfOurNetwork <- 20
SubNetworkSize <- 5
net<-BuildSemiCompleteNetwork(SizeOfOurNetwork,SubNetworkSize)
DM <- distances(net,v=V(net),to=V(net))  # distance mx
net<-SpatiallyCorrelatedLabels2(net,SubNetworkSize,2,1)
V(net)$MyAttribute

plot(net, edge.arrow.size=.5,  vertex.label.dist=1.5, vertex.size=5+2*V(net)$MyAttribute,main="original attributes")

#LISAnet <- MyLISA(net)
#plot(LISAnet, edge.arrow.size=.5, vertex.label.color="black", vertex.label.dist=1.5,
#vertex.size=0.25*(0.01-min(V(LISAnet)$LISAstat)+V(LISAnet)$LISAstat),main="LISA values")

# permute the network and re-evaluate
i <- 1
V(net)$PermedAttribute <- FixedPerm(V(net)$MyAttribute,i)
V(net)$PermedAttribute
#plot(net, edge.arrow.size=.5,  vertex.label.dist=1.5, vertex.size=5+2*V(net)$PermedAttribute,main=paste("permuted attributes with node",i," fixed"))
#LISAnet2 <- MyLISA(net)

net2 <- CopyNet(net)
plot(net2, edge.arrow.size=.5,  vertex.label.dist=1.5, vertex.size=5+2*V(net2)$MyAttribute,main="attributes of permuted network")

NPerms <- 20
PValue <- rep(0,SizeOfOurNetwork)
LStat <- rep(0,NPerms)

f <- function(j){
    V(net2)$MyAttribute <- FixedPerm(V(net2)$MyAttribute,index)
    LStat[j] <- LisaNetForNodeI(net2,index)
}
for (index in 1:SizeOfOurNetwork)
{
  cat("\n",index)
  # Evaluate Lisa stat for node i
  LStat=mclapply(1:NPerms, f)
  LStat=unlist(LStat)
  #boxplot(LStat)
  abline(h=V(LISAnet)$LISAstat[index],col="red")
  PValue[index] <- sum(LStat>V(LISAnet)$LISAstat[index])/NPerms
}
(PValue)

#boxplot(PValue[1:25])
#boxplot(PValue[26:200])
OverExp <- c(rep(0,25),rep(1,175))
PlotData <- as.data.frame(cbind(PValue,OverExp))

ggplot(PlotData, aes(x = factor(OverExp), y = PValue)) +
  geom_violin() +
   geom_point() +
   xlab("Group") +
  ylab("Permutation test P-value") +
   theme_grey(base_size = 24, base_family = "Times")
```

Let's repeat that, but on more realistic graphs, This time we recording the lowest p-value across all nodes for each (permuted) graph.
Let's base the analysis on the alzheimer's network from Bryan, but using subgraphs chosen by chosing a node at random and then iteratively adding neighbors until we have 100 loci

```{r AlzBasedTests, eval=FALSE}
FullNetwork <- alz
plot(FullNetwork,vertex.label.cex=0.2,edge.arrow.size=0.3,vertex.size=2)

set.seed(23)
# work with a subgraph of 50 nodes
SubGraphSize <- 50

# generate the subgraph 
# Need a new version of this function that passes the focal node as a parameter
NewSub <- BuildSubGraph(FullNetwork,SubGraphSize)

plot(NewSub,vertex.label=V(NewSub)$MyID,vertex.label.cex=0.8,edge.arrow.size=0.3)
cat("\n#components: ",count_components(NewSub))
cat("\nFocalNode: ",FocalNode)

newnet <- SpatiallyCorrelatedLabels(NewSub,FocalNode)
plot(newnet, vertex.color=ifelse(V(newnet)==FocalNode,"blue","red"), vertex.label.dist=1.5,vertex.size=4+10*V(newnet)$MyAttribute)
  
cat("\nMorans-I: ",MoransI(distances(newnet,v=V(newnet),to=V(newnet)),V(newnet)$MyAttribute),"    expectation= ",-1/(length(V(newnet)$MyAttribute)-1))
  
# Does it look significant?
#Results2 <-PermutationTest(newnet,10,1)
#hist(Results2,breaks=20,col="grey")
#abline(v=MoransI(DM,V(newnet)$MyAttribute),col="red")
#cat("\nObserved Mortan-s-I value is ",MoransI(DM,V(newnet)$MyAttribute))
```

Here's a version of BuildSubGraph that takes the focal node as an argument

```{r FocalSubGraph}

BuildFocalSubGraph <- function(g, subgsize, FocalNodeID) {
  size <- 0
  # label the vertices with our own IDs, that can be copied over
  V(g)$MyID <- 1:vcount(g)
  
  # pick a random vertex to start with (this may go wrong if the graph isn't fully-connected)
  # subg <- sample(1:vcount(g), 1)  # this was how the old function did it
  subg <- FocalNodeID
  size <- length(subg)
  itcount<-1
  #cat("\nIteration: ",itcount,"  Subgraph size=",size,"  vertices: ",subg)
  while  ((size < subgsize)&&(itcount<1e6))
  {
    itcount <- itcount+1
    # pick a random vertex in the graph currently
    if (length(subg)==1){
      v <- subg
    }
    else
    {
      v <- sample(subg,1)
    }
    # get the neighbors of v
    v1 <- neighbors(g, V(g)[v], mode="all")
    # sample one of those
    if (length(v1)>0){
      if (length(v1)==1){
        v2 <- v1
      }else{
         v2 <- sample(v1, 1)
      }
      # add it to subg
      if (!(v2 %in% subg)){
        #cat("\nFocal vertex: ",v,"  nbrs: ",v1,"  sampled: ",v2)
        subg <- c(subg, v2)
      } 
      size <- length(subg)
    }
    #cat("\nIteration: ",itcount,"  Subgraph size=",size,"  vertices: ",subg)
    
    # debugging:
    compcount <- count_components(subgraph(g,vids=subg))
    if (compcount>1){
      cat("\n***Error. Disconnected subgraph")
      tempsub <- subgraph(g,vids=subg)
      #plot(tempsub,vertex.label=V(tempsub)$MyID,vertex.label.cex=0.8,edge.arrow.size=0.3)
     return(tempsub)
    }
  }
  Sub <- subgraph(g,vids=subg)
  if (size==subgsize){
    return(Sub)
  }else{
    cat("\nSub-graph routine failed: too many iterations")
    return(NA)
  }
}   
```


Repeating the above on alz-based graphs using the new subgraph function.
We record the lowest p-value across all nodes for each (permuted) graph.
We use subgraphs chosen by chosing a node at random and then iteratively adding neighbors until we have 100 loci

```{r AlzBasedTests2}
#FullNetwork <- alz
FullNetwork <- S

SpatiallyCorrelatedLabels3 <- function(ntwk,CentralNode,NeighborhoodDistance){
  # label CentralNode and all its near neighbors 5, and everything else Normal(0,1).

  # Find the focal node
  CenterNode <- which(V(ntwk)$MyID==CentralNode)
  
  # The focal node will always be the first node - check to make sure...
  #if (V(ntwk)$MyID[1] != CentralNode){
  #  # throw an error
  #  cat("\nCentral node was incorrect in SpatiallyCorrelatedLabels3()")
  #  cat("\nCentralNode")
  #  break;
  #}else{
  #  CentralNode <- 1
  #}
  
  NbrCount<- -1  # Don't want to count the Focal Node itself
  for (i in 1:length(V(ntwk))){
    if (distances(ntwk,CenterNode,i) <= NeighborhoodDistance){
      V(ntwk)$MyAttribute[i] <- 5
      NbrCount <- NbrCount +1
    }
    else
      V(ntwk)$MyAttribute[i] <- rnorm(1,0,1)
  }
  cat("\nThere are ",NbrCount, " nodes within  distance of ",NeighborhoodDistance)
  return (ntwk)
}

set.seed(23456)
# work with a subgraph of 50 nodes
SubGraphSize <- 10

# generate the subgraph 
Focus <- sample(1:vcount(FullNetwork), 1)  # this was how the old function did it

NewSub <- BuildFocalSubGraph(FullNetwork,SubGraphSize,Focus)

plot(NewSub,vertex.label=V(NewSub)$MyID,vertex.label.cex=0.8,edge.arrow.size=0.3,main="Subgraph")
cat("\nSubgraph #components: ",count_components(NewSub))
cat("\nFocalNode: ",Focus)

newnet <- SpatiallyCorrelatedLabels3(NewSub,Focus,2)
plot(newnet, vertex.color=ifelse(V(newnet)==Focus,"blue","red"), vertex.label.dist=2.5,edge.arrow.size=0.3,vertex.label.cex=ifelse(V(newnet)==Focus,1,0.4),vertex.size=ifelse(V(newnet)==Focus,10,5))
  
cat("\nMorans-I: ",MoransI(distances(newnet,v=V(newnet),to=V(newnet)),V(newnet)$MyAttribute),"    expectation= ",-1/(length(V(newnet)$MyAttribute)-1))


# Generate some attribute values - use FocalValue for the focal node and it's near neighbors; Unif(0,1) everywhere else


newnet <- SpatiallyCorrelatedLabels3(newnet,Focus,2)

plot(newnet, vertex.color=ifelse(V(newnet)==Focus,"blue","red"), vertex.label.dist=2.5,edge.arrow.size=0.3,vertex.label.cex=0.2,vertex.size= as.numeric(ifelse(V(newnet)$MyAttribute < 0,"0.1",2+V(newnet)$MyAttribute)),main="Subgraph showing attribute values")

  
# Do a permutation test to see if the results when testing using Morans I look significant?
# The nextline causes an error - to be fixed
#Results3 <-PermutationTest(newnet,100,1)
#hist(Results3,breaks=20,col="grey")
#abline(v=MoransI(DM,V(newnet)$MyAttribute),col="red")
#cat("\nObserved Moran-s-I value is ",MoransI(DM,V(newnet)$MyAttribute))
```



Now let's do some permutation tests on this graph using local Moran's I.
The permutation test is now written as a function.

```{r LocalMoran}


CopyNet2 <- function(ntwk){
  return (ntwk)
}

set.seed(45642)  # Error when seed is 45644 - fixed
NPerms <- 250

# work with a subgraph 
SubGraphSize <- 250

# generate the subgraph 
FocalNode <- sample(1:vcount(FullNetwork), 1)  # this was how the old function did it

NewSub <- BuildFocalSubGraph(FullNetwork,SubGraphSize,FocalNode)
WhichIsFocal <- which(V(NewSub)$MyID==FocalNode)


# Generate some attribute values
newnet <- SpatiallyCorrelatedLabels3(NewSub,FocalNode,1)
NewSub <- CopyNet2(newnet) 


plot(NewSub,vertex.label=V(NewSub)$MyID,vertex.label.cex=0.4,edge.arrow.size=0.1,main="Subgraph")

# we are using an old  version of LisaNetForNodeI() 
# and re-using some of the code above

OldLisaNetForNodeI <- function(ntwk,IthNode){
  V(ntwk)$StandardizedAttribute <- rep(-9,length(V(ntwk)))
  AttributeMean <- mean(V(ntwk)$MyAttribute)
  # Standardize node labels - subtract the mean attribute value from all attributes
  V(ntwk)$StandardizedAttribute <- V(ntwk)$MyAttribute - AttributeMean
  for (i in IthNode:IthNode){  # This loop can be removed. It's a legacy of code I borrowed from a Moran's-I function
    L <- 0
    for (j in 1:length(V(ntwk))){
      if (distances(ntwk,v=V(ntwk)[i],to=V(ntwk)[j]) == 1){
        # they are neighbors - for now we only consider directly adjacent nodes
        L <- L + V(ntwk)$StandardizedAttribute[j] # there would also be a weight term here if we wanted a more general definition of 'neighbor'
      }
    }
    L <- L * V(ntwk)$StandardizedAttribute[i]
  }
  return (L)
}

plot(NewSub, edge.arrow.size=.2,  vertex.size= as.numeric(ifelse(V(NewSub)$MyAttribute < 0,"0.1",2+V(NewSub)$MyAttribute)),main="original attributes",vertex.label.cex=0.2)

# OBserved LStat value:
ObservedLISA <- OldLisaNetForNodeI(NewSub,WhichIsFocal)

# Function to permute the attribute values and recalculate Local Moran's I
#PValue <- rep(0,SizeOfOurNetwork)
LStat <- rep(0,NPerms)

PermuteAndRecalculate <- function(ntwk,Focus){
  WhichIsFocal <- which(V(ntwk)$MyID==Focus)
  V(ntwk)$PermedAttribute <- FixedPerm(V(ntwk)$MyAttribute,WhichIsFocal)
  V(ntwk)$MyAttribute <- V(ntwk)$PermedAttribute

  #plot(ntwk, edge.arrow.size=.2,  vertex.label.cex=0.2, vertex.size= as.numeric(ifelse(V(ntwk)$PermedAttribute < 0,"0.1",2+V(ntwk)$PermedAttribute)),main="attributes of sample permuted network")
  
  LMI <- OldLisaNetForNodeI(ntwk,WhichIsFocal)
  #cat("\nLMI:",LMI)
  return(LMI)
}

#(head(V(net2)$MyAttribute))
for (kk in 1:NPerms){
  LStat[kk] <- PermuteAndRecalculate(NewSub,FocalNode)
}
#(head(V(NewSub)$MyAttribute))
PValue <- sum(LStat>ObservedLISA)/NPerms
cat("\nNewPValue: ",PValue)
boxplot(LStat)
abline(h=ObservedLISA,col="red",main="New distribution of LISA Stat")
cat("\nLISA stat value for focal node: ",ObservedLISA)
NewLStat <- LStat



# The following is old code to check we still get the same answer now that we have re0written thingds...

# permute the network and re-evaluate
net2 <- CopyNet2(NewSub)
i <- FocalNode
WhichIsFocal <- which(V(net2)$MyID==FocalNode)
V(net2)$PermedAttribute <- FixedPerm(V(net2)$MyAttribute,WhichIsFocal)
plot(net2, edge.arrow.size=.2,  vertex.label.cex=0.2, vertex.size= as.numeric(ifelse(V(net2)$PermedAttribute < 0,"0.1",2+V(net2)$PermedAttribute)),main="attributes of sample permuted network")


ff <- function(j){
    V(net2)$MyAttribute <- FixedPerm(V(net2)$MyAttribute,WhichIsFocal)
    LStat[j] <- OldLisaNetForNodeI(net2,WhichIsFocal)
}

for (i in WhichIsFocal:WhichIsFocal)
{
  #cat("\n",index)
  # Evaluate Lisa stat for node i
  LStat=mclapply(1:NPerms, ff)
  LStat=unlist(LStat)
  PValue[WhichIsFocal] <- sum(LStat>ObservedLISA)/NPerms
}

boxplot(LStat)
abline(h=ObservedLISA,col="red",main="Distribution of LISA Stat")
cat("\nLISA stat value for focal node: ",ObservedLISA)


cat("\nKS test result:\n")
(ks.test(NewLStat,LStat))
ks <- ks.test(NewLStat,LStat)
plot(ecdf(LStat), 
     xlim = range(c(LStat, NewLStat)), 
     col = "blue",pch='.')
plot(ecdf(NewLStat), 
     add = TRUE, 
     lty = "dashed",
     col = "red",pch='.')
 

#OverExp <- c(rep(0,25),rep(1,175))
#PlotData <- as.data.frame(cbind(PValue,OverExp))
#ggplot(PlotData, aes(x = factor(OverExp), y = PValue)) +
#  geom_violin() +
#   geom_point() +
#   xlab("Group") +
#  ylab("Permutation test P-value") +
#   theme_grey(base_size = 24, base_family = "Times")


```

Now let's try to do repeated permutation tests and look at the distribution of the resulting p-values....
```{r FullPerm}

library(foreach)
set.seed(95)  #(45643)  
NPerms <- 100 #250
NRepeats <- 20
PVals <- rep(-9,NRepeats)
OrigLISAs <- rep(-9,NRepeats)
nNeighbors <- rep(0,NRepeats)
SubGraphSize <- 500 ## 500

SpatiallyCorrelatedLabels4 <- function(ntwk,CentralNode,NeighborhoodDistance,FocalMean,RemainingMean){
  # label CentralNode and all its near neighbors Normal(FocalMean,1), and everything else Normal(RemainingMean,1).

  # Find the focal node
  CenterNode <- which(V(ntwk)$MyID==CentralNode)
  
  NbrCount<- -1  # Don't want to count the Focal Node itself
  for (i in 1:length(V(ntwk))){
    if ((distances(ntwk,CenterNode,i) <= NeighborhoodDistance)&&(NbrCount<50000)){
      V(ntwk)$MyAttribute[i] <- rnorm(1,FocalMean,1)
      NbrCount <- NbrCount +1
    }
    else
      V(ntwk)$MyAttribute[i] <- rnorm(1,RemainingMean,1)
  }
  cat("\nThere are ",NbrCount, " nodes within  distance of ",NeighborhoodDistance)
  nNeighbors[mm] <<-NbrCount
  return (ntwk)
}

foreach (mm=1:NRepeats) %do% {
  
  # generate the subgraph 
  FocalNode <- sample(1:vcount(FullNetwork), 1)  # this was how the old function did it

  NewSub <- BuildFocalSubGraph(FullNetwork,SubGraphSize,FocalNode)
  WhichIsFocal <- which(V(NewSub)$MyID==FocalNode)
   
  # Generate some attribute values
  newnet <- SpatiallyCorrelatedLabels4(NewSub,FocalNode,2,3,2)
  NewSub <- CopyNet2(newnet) 
 
  #plot(NewSub, edge.arrow.size=.2,  vertex.size= as.numeric(ifelse(V(NewSub)$MyAttribute < 0,"0.1",2+V(NewSub)$MyAttribute)),main=c("original atts. Iter: ",mm),vertex.label.cex=0.2)

  ObservedLISA <- OldLisaNetForNodeI(NewSub,WhichIsFocal)
  OrigLISAs[mm] <- ObservedLISA
  LStat <- rep(0,NPerms)
  
  for (kk in 1:NPerms){
    LStat[kk] <- PermuteAndRecalculate(NewSub,FocalNode)
  }
  #(head(V(NewSub)$MyAttribute))
  PVals[mm] <- sum(LStat>ObservedLISA)/NPerms
  cat("\nNewPValue: ",PVals[mm],"   orig. LisaStat: ",ObservedLISA)
  if (1){ #(PVals[mm]>0){
      plot(NewSub, edge.arrow.size=.2,  vertex.size= as.numeric(ifelse(V(NewSub)$MyAttribute < 2,"2",2+3*V(NewSub)$MyAttribute)),main=c("original atts. Iter: ",mm," p=",PVals[mm]),vertex.label.cex=0.2)

  }

}

(boxplot(PVals))
plot(nNeighbors,PVals)


```

Let's look at how the subgraphs resemble the overall super network in terms of various summary statistics:
```{r,   NetworkSummaries}
y <- degree_distribution(FullNetwork)[1:150]
x <- 1:length(y)
lo <- loess(y~x)
plot(x,y,pch='.',main="Supernetwork Degree Distn (truncated at 150")
lines(predict(lo), col='red', lwd=2)


for (ii in 1:5){
  FocalNode <- sample(1:vcount(FullNetwork), 1)  # this was how the old function did it
  NewSub <- BuildFocalSubGraph(FullNetwork,SubGraphSize,FocalNode)
  y <- degree_distribution(NewSub)
  x <- 1:length(y)
  lo <- loess(y~x)
  plot(x,y,pch='.',main="Example Subgraph Degree Distn")
  lines(predict(lo), col='red', lwd=2)
}

# add other graph diagnostics
cat("\nDiam. of S: ",diameter(FullNetwork))
cat("\nDiam.of e.g. subgraph: ",diameter(NewSub))
cat("\n\nAv. Path Len of S: ",average.path.length(FullNetwork))
cat("\nAv. Path Len of e.g. subgraph: ",average.path.length(NewSub))

# q: how best to compare two graphs in terms of their overall structure?
# a: use the graph edit distance, which is the minimum number of edge additions, deletions, and substitutions needed to transform one graph into the other.

```



***
&nbsp;
***
&nbsp;
***
&nbsp;
## Other methods

There are a variety of other methods out there to detect clustering. There is no need to pursue these now, but they are worth keeping in mind for future reference.

***

### Old notes: Tango's disease clusters (Tango 1999)

This method is written to look for 'regions' that have too many diseased indivs. 

As such, it does not really fit into our paradigm, because we have no regions. 

We could force it to (possibly?) work if we defined each gene node to be a 'region' in which there was just one potential case. 

It extends the original Tango method, in that it allows the size of the locally clustered region to be a parameter that is searched over.


***

### Kulldorf's local scan statistic (Kulldorf 1997)

This method is a local cluster detection algorithm that is implemented in the R package SaTScan.

The basic idea is to model, under the null, the occurrence of 'cases' as, in the simplest case, a Poisson process, and then attempt to detect whether this model, which assumes independence everywhere, fits. 

If it does not fit, clustering is assumed to be occurring. 

It is then extended to other 'null' measures as well, to account for uneven distributions of potential cases (I think), but I don't think we would need those ideas here.

The method originally tried to detected circular clusters of 'alike' values,but has been extended to non-circular shapes. Ideas like this would presumably be sensible if we used local neighborhoods rather than 'shapes'.

An important characteristic of this method is that when the null hypothesisis rejected it indicates the specific area of the map that causes the rejection.

Again, though, it is a region-based method, so not ideal for networks


***

### Besag and Newell's method (Besag et al. 1991)

This is a local cluster detection method that assumes the landscape is divided into regions. 

The paper talks about why the assumptions of Moran's-I might be invalid even when disease risk is $p$ for everyone and all are independent. 

Essentially, if the sizes of the regions varies, as it will, then the variance of the incidence rates will also vary from region to region. 

This actually induces spatial correlation in the disease rates even when the disease rate is constant from individual to individual. This won't be a problem in our context, because there is one potential 'case' in each region, so all 'disease' rates have the same distribution.


***

### Arbitrarily shaped multiple spatial cluster detection for case event data (Dematte?? et al. 2007)

This paper proposes a method for spatial cluster detection of case event data. 

Here, each case has a unique location in space.

Essentially, they look at a path that includes (just) all cases, and look at whether the average distance between points in this path is unusually low (suggesting the existence of clusters.)

**These ideas look like they might, with some thought, generalize to our setting. The network would just impose constraints on how the path is constructed.**

The paper also contains a nice (brief) overview of methods, referencing several recently proposed methods.


***