---
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 PMRplus is used to to obtain rank-based tests of accuracy. We also need the R package ``tidyverse'' to massage the data into the appropriate format for hypothesis testing. The references are given below.

Before the tests are run, we subtract the true value of Centrality from each of the centrality values for each of the 100 samples, then divide the difference by the true centrality, thus obtaining a percentage error. This operation gets all centrality measures on the same scale and allows comparisons with the true value.

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 PMCMRplus package are NPM_test ordered pairs test. The data need to be in long format and balanced across the groups.

The simulated data are in "/Users/monniemcgee/Dropbox/2019 Spring/FlowThroughCentrality"

```{r loadPackages}
#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 12% perturbation
```{r Nodes50deg6pert12}
setwd("/Users/monniemcgee/Dropbox/2019 Spring/FlowThroughCentrality/Simulations/Nodes50Deg6Pert12")
# 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_Pert12_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)
dataBW$Sample <- 1:100
dataBWLong <- gather(dataBW,key="CentLevel",value="Centrality",starts_with("Cent"), factor_key=TRUE) %>% select(Sample,Network,Measure,Centrality,CentLevel)

# Import closeness data and convert to long format
dataCL <- read.csv(file="Nodes50_AvgDeg06_UW_01_Pert12_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)
dataCL$Sample <- 1:100
dataCLLong <- gather(dataCL,key="CentLevel",value="Centrality",starts_with("Cent"), factor_key=TRUE) %>% select(Sample,Network,Measure,Centrality,CentLevel)

# Import flow between data and convert to long format
dataFB <- read.csv(file="Nodes50_AvgDeg06_UW_01_Pert12_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)
dataFB$Sample <- 1:100
dataFBLong <- gather(dataFB,key="CentLevel",value="Centrality",starts_with("Cent"), factor_key=TRUE) %>% select(Sample,Network,Measure,Centrality,CentLevel)

# Import Flowthrough data and convert to long format
dataFT <- read.csv(file="Nodes50_AvgDeg06_UW_01_Pert12_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)
dataFT$Sample <- 1:100
dataFTLong <- gather(dataFT,key="CentLevel",value="Centrality",starts_with("Cent"), factor_key=TRUE) %>% select(Sample,Network,Measure,Centrality,CentLevel)

# Import Stable betweenness data and convert to long format
dataSB <- read.csv(file="Nodes50_AvgDeg06_UW_01_Pert12_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)
dataSB$Sample <- 1:100
dataSBLong <- gather(dataSB,key="CentLevel",value="Centrality",starts_with("Cent"), factor_key=TRUE) %>% select(Sample,Network,Measure,Centrality,CentLevel)
```

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

## Nashimoto-Wright NPM-test for ordered means of non-normal data
## Select only one Level at a time - NPM test has no blocking variables
## We will do this for the first 5 levels.
n50deg6pert12lev1 <- subset(n50deg6pert12, CentLevel == "Cent1", select = c("Sample","Network","Measure","Centrality"))
npmLev1 <- NPMTest(Centrality~Measure,data=n50deg6pert12lev1)
n50deg6pert12lev2 <- subset(n50deg6pert12, CentLevel == "Cent2", select = c("Sample","Network","Measure","Centrality"))
npmLev2 <- NPMTest(Centrality~Measure,data=n50deg6pert12lev2)
n50deg6pert12lev3 <- subset(n50deg6pert12, CentLevel == "Cent3", select = c("Sample","Network","Measure","Centrality"))
npmLev3 <- NPMTest(Centrality~Measure,data=n50deg6pert12lev3)
n50deg6pert12lev4 <- subset(n50deg6pert12, CentLevel == "Cent4", select = c("Sample","Network","Measure","Centrality"))
npmLev4 <- NPMTest(Centrality~Measure,data=n50deg6pert12lev4)
n50deg6pert12lev5 <- subset(n50deg6pert12, CentLevel == "Cent5", select = c("Sample","Network","Measure","Centrality"))
npmLev5 <- NPMTest(Centrality~Measure,data=n50deg6pert12lev5)
npmLev1

npmLev2

npmLev3

npmLev4

npmLev5
```

```{r Graphics}
# Plotting only 5 levels due to overcrowding of plot.

n50deg6pert12all <- bind_rows(n50deg6pert12lev1,n50deg6pert12lev2,n50deg6pert12lev3,n50deg6pert12lev4,n50deg6pert12lev5, .id="CentLevel")
pdf(file="n50deg6pert12Lines.pdf")
ggplot(n50deg6pert12all, aes(x=Sample,y=Centrality)) + geom_line(aes(linetype=factor(Measure))) + facet_grid(CentLevel~.) + theme_bw() + scale_linetype_discrete(name="Centrality\nMeasure")
dev.off()

sumData <- n50deg6pert12all %>% group_by(CentLevel,Measure) %>% summarise(Mean=mean(Centrality),SD=sd(Centrality),Median=median(Centrality),MAD=mad(Centrality))
pdf(file="n50deg6pert12Summary.pdf")
ggplot(sumData,aes(x=CentLevel,y=Mean)) + geom_point(aes(shape=factor(Measure))) + theme_bw() + scale_shape_discrete(name="Centrality\nMeasure")
dev.off()

xtable(sumData)
```
