---
title: "FlowThrough2019"
author: "Monnie McGee"
date: "4/1/2019"
output: html_document
---

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

This document contains code used for the Flow-Through Centrality project, to test whether which method of computing centrality is the most stable. The R package coin is used to to obtain permutation-based tests of equality of variance. We also need the R package ``tidyverse'' to massage the data into the appropriate format for hypothesis testing. The references are given below.

Torsten Hothorn, Kurt Hornik, Mark A. van de Wiel, Achim Zeileis (2008). Implementing a Class of Permutation
  Tests: The coin Package. Journal of Statistical Software 28(8), 1-23. URL http://www.jstatsoft.org/v28/i08/

The functions we use from the coin package are fligner_test for the main hypothesis test. The data need to be in long format and balanced across the groups.

The simulated data are in "~/Dropbox/2019 Spring/FlowThroughCentrality/Simulations"


```{r loadPackages}
setwd("~/Dropbox/2019 Spring/FlowThroughCentrality")
#install.packages("coin") - use if you have not installed the package.
library(coin) # load the package in R 
library(tidyverse) # for tidying the data
library(xtable) # for obtaining tables in LaTex format
library(PMCMRplus) # For Nashimoto-Wright NPM-test for ordered means of non-normal datajo
```

# 50 Nodes with Six degrees and 3% perturbation
```{r Nodes50deg6pert3}
setwd("./Simulations/Nodes50Deg6Pert3")
# Create variable names for all data sets
varNames <- c("Network",paste0(c("Cent","Tier"), rep(1:50,each=2)))

# Import betweenness data and convert to long format
dataBW <- read.csv(file="Nodes50_AvgDeg06_UW_01_Pert03_Betweenness.csv",skip=4,nrows=101,header=FALSE)
colnames(dataBW) <- varNames
dataBW$Network<-str_replace_all(dataBW$Network, c(" " = "" ))
dataBW <- dataBW %>% mutate_if(is.numeric, .funs =function(x) 100*(x-x[1])/x[1])
dataBW <- dataBW[-1,]
dataBW$Measure <- rep("BW",100)
dataBWLong <- gather(dataBW,key="CentLevel",value="Centrality",starts_with("Cent"), factor_key=TRUE) %>% select(Network,Measure,Centrality,CentLevel)

# Import closeness data and convert to long format
dataCL <- read.csv(file="Nodes50_AvgDeg06_UW_01_Pert03_Closeness.csv",skip=4,nrows=101,header=FALSE)
colnames(dataCL) <- varNames
dataCL$Network<-str_replace_all(dataCL$Network, c(" " = "" ))
dataCL <- dataCL %>% mutate_if(is.numeric, .funs =function(x) 100*(x-x[1])/x[1])
dataCL <- dataCL[-1,]
dataCL$Measure <- rep("CL",100)
dataCLLong <- gather(dataCL,key="CentLevel",value="Centrality",starts_with("Cent"), factor_key=TRUE) %>% select(Network,Measure,Centrality,CentLevel)

# Import flow between data and convert to long format
dataFB <- read.csv(file="Nodes50_AvgDeg06_UW_01_Pert03_FlowBetween.csv",skip=4,nrows=101,header=FALSE)
colnames(dataFB) <- varNames
dataFB$Network<-str_replace_all(dataFB$Network, c(" " = "" ))
dataFB <- dataFB %>% mutate_if(is.numeric, .funs =function(x) 100*(x-x[1])/x[1])
dataFB <- dataFB[-1,]
dataFB$Measure <- rep("FB",100)
dataFBLong <- gather(dataFB,key="CentLevel",value="Centrality",starts_with("Cent"), factor_key=TRUE) %>% select(Network,Measure,Centrality,CentLevel)

# Import Flowthrough data and convert to long format
dataFT <- read.csv(file="Nodes50_AvgDeg06_UW_01_Pert03_Flowthrough.csv",skip=4,nrows=101,header=FALSE)
colnames(dataFT) <- varNames
dataFT$Network<-str_replace_all(dataFT$Network, c(" " = "" ))
dataFT <- dataFT %>% mutate_if(is.numeric, .funs =function(x) 100*(x-x[1])/x[1])
dataFT <- dataFT[-1,]
dataFT$Measure <- rep("FT",100)
dataFTLong <- gather(dataFT,key="CentLevel",value="Centrality",starts_with("Cent"), factor_key=TRUE) %>% select(Network,Measure,Centrality,CentLevel)

# Import Stable betweenness data and convert to long format
dataSB <- read.csv(file="Nodes50_AvgDeg06_UW_01_Pert03_StableBetween.csv",skip=4,nrows=101,header=FALSE)
colnames(dataSB) <- varNames
dataSB$Network<-str_replace_all(dataSB$Network, c(" " = "" ))
dataSB <- dataSB %>% mutate_if(is.numeric, .funs =function(x) 100*(x-x[1])/x[1])
dataSB <- dataSB[-1,]
dataSB$Measure <- rep("SB",100)
dataSBLong <- gather(dataSB,key="CentLevel",value="Centrality",starts_with("Cent"), factor_key=TRUE) %>% select(Network,Measure,Centrality,CentLevel)

## Bind all data sets together in long format
n50deg6pert3 <- bind_rows(dataFTLong, dataFBLong, dataCLLong, dataBWLong, dataSBLong)
n50deg6pert3$Measure <- ordered(n50deg6pert3$Measure, levels=c("FT","FB","CL","BW","SB"))

## Jonckheere-Terpstra test for ordered alternatives
## Select only one Level at a time - JT test has no blocking variables
## We will do this for the first 5 levels.
## An example: ed_exp4 <- subset(education, Region == 2, select = c("State","Minor.Population","Education.Expenditures"))
n50deg6pert3lev1 <- subset(n50deg6pert3, CentLevel == "Cent1", select = c("Network","Measure","Centrality"))
npmLev1 <- NPMTest(Centrality~Measure,data=n50deg6pert3lev1)

# Coefficient of variation function (can't believe I had to write this!)  
cv <- function(x){
  cv=100*(sd(x)/mean(x))
}

# There are 100 Levels, therefore, this plot is too messy. Need to think of something else.
# pdf(file="n50deg6pert3Lines.pdf")
# ggplot(n50deg6pert3,aes(x=Network, y=Centrality, colour=Measure)) + geom_line() + facet_grid(Level~.,scales="free")
# dev.off()


sumData <- n50deg6pert3 %>% group_by(CentLevel,Measure) %>% summarise(Mean=mean(Centrality),SD=sd(Centrality),CV=cv(Centrality))
ggplot(sumData,aes(x=CentLevel,y=CV,colour=Measure)) + geom_line()

# xtable(sumData) 

```

The Fligner-Kileen test is a test for homogeneity of variance in multiple independent samples. Here, we treat each measure of centrality as a different sample, and test to determine which measure is the most stable; in other words, has the lowest variance. 

The modified FK test jointly ranks the absolute differences of each observation from the median observation for that group (or centrality measure, in this case). It assigns increasing scores a quantile of the Normal distribution function, based on the ranks of all observations. Its test statistic has an asymptotic chi-square distribution with k-1 degrees of freedom, where k is the number of independent samples.

```{r HypothesisTesting}
# Need to do this test by Levels
# Multiple testing correction?
fligner_test(Centrality~Measure|Level,data=tier6n50deg6pert3)

sumData <- tier6n50deg6pert3 %>% group_by(Level,Measure) %>% summarise(Mean=mean(Centrality), SD=sd(Centrality),CV=cv(Centrality),Median=median(Centrality),MAD=mad(Centrality))

## Approximative (Monte Carlo) Lepage test
## Hollander and Wolfe (1999, p. 172)
# lepage_trafo <- function(y) cbind("Location" = rank_trafo(y), "Scale" = ansari_trafo(y))

## Why was the null hypothesis rejected?
## Note: maximum statistic instead of quadratic form
ltm <- independence_test(Centrality ~ Measure, data = tier6n50deg6pert3,
                  distribution = approximate(B = 10000),
                  ytrafo = function(data)
                      trafo(data, numeric_trafo = fligner_trafo))

pvalue(ltm, method="step-down")

# Plot the FKscores
tier6n50deg6pert3$FK <- fligner_trafo(tier6n50deg6pert3$Centrality)
ggplot(tier6n50deg6pert3,aes(x=Sample, y=FK, colour=Measure)) + geom_line() + facet_grid(CentLevel~.,scales="free")

```

