---
title: "Analysis of the Stability of Flow Through Centrality"
author: "Monnie McGee"
date: "2/14/2018"
output: html_document
bibliography: flowthrough.bib
---

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

## R Markdown

This document shows an analysis of preliminary data from a comparison of flow through centrality with betweenness centrality measures for three well-known networks. The networks are Zachary's Karate Club, Dolphin Families, and 9-11 Hijackers.


```{r SetWD, echo=FALSE}
# All files are in the following working directory and subdirectories.
# setwd("/Users/Monnie/Dropbox/2018 Spring Work/FlowThroughCentrality")
# Karate Club Data Subdirectory: "First_Pilot_KarateClub"                         
# 9-11 Hijackers Subdirectory: "First_Pilot_Test_9-11_Hijackers""       
# Dolphins Subdirectory: "First_Pilot_Test_Dolphins" 
require(openxlsx)
# Reading and Cleaning Data for Karate Club Betweenness Data at 3%
karate3bt <- read.xlsx("First_Pilot_KarateClub/KarateClub.betweenness_03.xlsx",startRow=3,rows=3:37)
varnames.bt <- c("Node","Bet.Orig","Tier.Orig","AvgBT","MaxBT","MinBT","AvgTier","HighTier","LowTier","BT1","Tier1","BT2","Tier2","BT3","Tier3","BT4","Tier4","BT5","Tier5","BT6","Tier6","BT7","Tier7","BT8","Tier8","BT9","Tier9","BT10","Tier10")
names(karate3bt) <- varnames.bt
karate3bt.cen <- karate3bt[,c(1,2,seq(10,28,by=2))]
karate3bt.tier <- karate3bt[,c(1,3,seq(11,29,by=2))]

# For karate club flow data at 3%
karate3fl <- read.xlsx("First_Pilot_KarateClub/KarateClub.flow_through_03.xlsx",startRow=3,rows=3:37)
varnames.fl <- c("Node","Flow.Orig","Tier.Orig","AvgFT","MaxFT","MinFT","AvgTier","HighTier","LowTier","FT1","Tier1","FT2","Tier2","FT3","Tier3","FT4","Tier4","FT5","Tier5","FT6","Tier6","FT7","Tier7","FT8","Tier8","FT9","Tier9","FT10","Tier10")
names(karate3fl) <- varnames.fl
karate3fl.cen <- karate3fl[,c(1,2,seq(10,28,by=2))]
karate3fl.tier <- karate3fl[,c(1,3,seq(11,29,by=2))]
```

First we divide each column pertaining to betweenness centrality by the maximum value for that column in order to facilitate comparisons between BC and FT, which is on a scale of 0 to 1.

```{r Scaling}
# Scale BC karate data by maximum
karate3bt.cen.scaled <- apply(karate3bt.cen[,2:12],2,function(x){x/max(x)})
karate3bt.cen.scaled <- as.data.frame(karate3bt.cen.scaled)
karate3bt.cen.scaled['Node'] <- karate3bt.cen$Node # add Node column
```

Each of the three networks was perturbed at 3%, 6%, 9% and 12%. For each of the perturbations, betweenness and flow through centralities were calculated 10 times. The plots below are line plots of the centralities for all 10 runs from the 3% perturbation of the Karate Club Data. The x-axis is the node. The vertical axis shows scaled betweenness centralities (first plot) in solid lines and flowthrough centralities (second plot) in dashed lines.

```{r BetweenPlot}
require(ggplot2)
require(reshape)
df.bet <- melt(karate3bt.cen.scaled,  id.vars = 'Node', variable.name = 'Value')
ggplot(df.bet, aes(Node,value)) + ylim(0,1) + geom_line(linetype=1,aes(y = value, colour = variable))
df.tier <- melt(karate3bt.tier, id.vars='Node', variable.name = "Value")
ggplot(df.tier, aes(Node,value)) + geom_line(linetype=1, aes(y=value, colour = variable))

df.fl <- melt(karate3fl.cen,  id.vars = 'Node', variable.name = 'Value')
ggplot(df.fl, aes(Node,value)) + ylim(0,1) + geom_line(linetype=2, aes(y=value, colour = variable))
df.tier.fl <- melt(karate3fl.tier, id.vars='Node', variable.name = "Value")
ggplot(df.tier.fl, aes(Node,value)) + geom_line(linetype=2, aes(y=value, colour = variable))
```

We wish to test the hypothesis that flow through centrality is more stable than betweenness centrailty. This amounts to testing the null hypothesis that the variability in flow through centrality is greater than or equal to the variability in betweenness centrality versus the alternative hypothesis that the varability in flow through centrality is less than that of betweenness centrality. However, we also have to keep in mind that centrality is measured by node, and that the nodes are the same for the networks we are measuring, no matter what method we use to measure the importance of a node. In other words, the nodes are paired. Statistically, this means that centrality measures for the same network should be correlated.

A basic statistical test used to examine differences in paired data is the paired t-test. However, the paired t-test is a test for differences in the mean. Here, we are interested in differences in variability between paired observations. Information on tests for differences in variability when the samples are independent is found in most introductory statistics textbooks. Information on tests for differences in variability when the samples are correlated is a relatively new research topic [ @seier; @bonett; @wilcox].

For the purposes of this document, we will use methods developed by @wilcox, since they are readily available in the R package BLAH (need to insert this information). 

## Note to self: insert detailed information about the method used here.

```{r RMD Rep1}
### RMD ###
# Null hypothesis is ratio of flow to betweeness <= 1.
# Alternative hypothesis is ratio of flow to betweenness is > 1.
dev.flow <- karate3fl.cen$FT1 - median(karate3fl.cen$Flow.Orig)
dev.bet <- karate3bt.cen.scaled$BT1 - median(karate3bt.cen.scaled$Bet.Orig)
R <- 9999
all <- c(dev.flow,dev.bet)
k <- 1:length(all)
reps <- numeric(R)
RMD <- (sum(abs(dev.flow))/length(dev.flow))/(sum(abs(dev.bet))/length(dev.bet))
RMD
for (i in 1:R) {
m <- sample(k, size=length(dev.flow), replace=FALSE)
p.flow <- all[m]
p.bet <- all[-m]
reps[i] <- (sum(abs(p.flow))/length(p.flow))/(sum(abs(p.bet))/length(p.bet))
}
# calculate one-tailed p-value
if (RMD > 1) pvalue <- mean(c(RMD , reps) >= RMD) + mean(c(RMD , reps) <= 1/RMD)
# if (RMD < 1) pvalue <- mean(c(RMD , reps) <= RMD) + mean(c(RMD , reps) >= 1/RMD)
# if (RMD == 1) pvalue <- 1
pvalue
```

Create the RMD function for other columns
```{r CalcRMD, echo=FALSE}
RMD <- function(x,y) {
  #x and y are vectors
 dev.x <- x - median(x)
 dev.y <- y - median(y)
 RMD <- (sum(abs(dev.x))/length(dev.x))/(sum(abs(dev.y))/length(dev.y))
 RMD
}

```

```{r RMD Rep2, echo=FALSE}
### RMD ###
# Null hypothesis is ratio of flow to betweeness <= 1.
# Alternative hypothesis is ratio of flow to betweenness is > 1.
dev.flow <- karate3fl.cen$FT2 - median(karate3fl.cen$FT2)
dev.bet <- karate3bt.cen.scaled$BT2 - median(karate3bt.cen.scaled$BT2)
R <- 9999
all <- c(dev.flow,dev.bet)
k <- 1:length(all)
reps <- numeric(R)
RMD <- (sum(abs(dev.flow))/length(dev.flow))/(sum(abs(dev.bet))/length(dev.bet))
RMD
for (i in 1:R) {
m <- sample(k, size=length(dev.flow), replace=FALSE)
p.flow <- all[m]
p.bet <- all[-m]
reps[i] <- (sum(abs(p.flow))/length(p.flow))/(sum(abs(p.bet))/length(p.bet))
}
# calculate two-tailed p-value
if (RMD > 1) pvalue <- mean(c(RMD , reps) >= RMD) + mean(c(RMD , reps) <= 1/RMD)
# if (RMD < 1) pvalue <- mean(c(RMD , reps) <= RMD) + mean(c(RMD , reps) >= 1/RMD)
# if (RMD == 1) pvalue <- 1
pvalue
```

# References

