---
title: "Hijacker Data Analysis for 9% Perturbations"
author: "Monnie McGee"
date: "July 12, 2022"
output: html_document
---

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

This document contains code used to analyse real data from the 9-11 Hijacker network for the Flow-Through Centrality project. We want to test whether which method of computing centrality is the most stable as applied to real data. 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 observed centrality values for each of the 10 samples, multiply by 2, then divide by the difference of the true centrality and the observed 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.

.csv files associated with this markdown document are in "~/Code/HijackerData". Original data were given in Excel files, which are in the directory "~/Code/HijackerData/ExcelFiles".

```{r loadPackages}
# setwd("./Code/HijackerData")
#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 data
```

RMarkdown file to enter Hijacker Data with 9% perturbation and calculate relative percentage error for each centrality measure. 

# Create variable names for all data sets
```{r VarNames}
varNames <- c("Network",paste0(c("Cent","Tier"), rep(1:62,each=2)))
```
# Import betweenness data and convert to long format
```{r ImportBetwenness}
dataBW <- read.csv(file="Hijack_9percent_betweenness.csv",skip=5,nrows=11,header=FALSE)
colnames(dataBW) <- varNames
dataBW$Network<-str_replace_all(dataBW$Network, c(" " = "" ))
test <- dataBW %>% mutate_at(vars(matches("Cent")), .funs = function(x) is.infinite(x))
length(test[test == TRUE])
dataBW <- dataBW %>% mutate_at(vars(matches("Cent")), .funs =function(x) (2*(x-x[1]))/(abs(x)+abs(x[1])))
dataBW <- dataBW %>% mutate_at(vars(matches("Cent")), .funs = function(x) ifelse(is.nan(x), 0, x))
dataBW <- dataBW[-1,] # delete first row of true values
dataBW$Measure <- rep("BW",10)
dataBW$Sample <- 1:10
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
```{r ImportCloseness}
dataCL <- read.csv(file="Hijack_9percent_Closeness.csv",skip=5,nrows=11,header=FALSE)
colnames(dataCL) <- varNames
dataCL$Network<-str_replace_all(dataCL$Network, c(" " = "" ))
test <- dataCL %>% mutate_at(vars(matches("Cent")), .funs = function(x) is.infinite(x))
length(test[test == TRUE])
dataCL <- dataCL %>% mutate_at(vars(matches("Cent")), .funs =function(x) (2*(x-x[1]))/(abs(x)+abs(x[1])))
dataCL <- dataCL %>% mutate_at(vars(matches("Cent")), .funs = function(x) ifelse(is.nan(x), 0, x))
dataCL <- dataCL[-1,] # delete first row of true values
dataCL$Measure <- rep("CL",10)
dataCL$Sample <- 1:10
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
```{r ImportFlowBetween}
dataFB <- read.csv(file="Hijack_9percent_FlowBetween.csv",skip=5,nrows=11,header=FALSE)
colnames(dataFB) <- varNames
dataFB$Network<-str_replace_all(dataFB$Network, c(" " = "" ))
test <- dataFB %>% mutate_at(vars(matches("Cent")), .funs = function(x) is.infinite(x))
length(test[test == TRUE])
dataFB <- dataFB %>% mutate_at(vars(matches("Cent")), .funs =function(x) (2*(x-x[1]))/(abs(x)+abs(x[1])))
dataFB <- dataFB %>% mutate_at(vars(matches("Cent")), .funs = function(x) ifelse(is.nan(x), 0, x))
dataFB <- dataFB[-1,] # delete first row of true values
dataFB$Measure <- rep("FB",10)
dataFB$Sample <- 1:10
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
```{r ImportFlowthrough}
dataFT <- read.csv(file="Hijack_9percent_Flowthrough.csv",skip=5,nrows=11,header=FALSE)
colnames(dataFT) <- varNames
dataFT$Network<-str_replace_all(dataFT$Network, c(" " = "" ))
test <- dataFT %>% mutate_at(vars(matches("Cent")), .funs = function(x) is.infinite(x))
length(test[test == TRUE])
dataFT <- dataFT %>% mutate_at(vars(matches("Cent")), .funs =function(x) (2*(x-x[1]))/(abs(x)+abs(x[1])))
dataFT <- dataFT %>% mutate_at(vars(matches("Cent")), .funs = function(x) ifelse(is.nan(x), 0, x))
dataFT <- dataFT[-1,] # delete first row of true values
dataFT$Measure <- rep("FT",10)
dataFT$Sample <- 1:10
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
```{r ImportStableBetweenness}
dataSB <- read.csv(file="Hijack_9percent_StableBetween.csv",skip=5,nrows=11,header=FALSE)
colnames(dataSB) <- varNames
dataSB$Network<-str_replace_all(dataSB$Network, c(" " = "" ))
test <- dataSB %>% mutate_at(vars(matches("Cent")), .funs = function(x) is.infinite(x))
length(test[test == TRUE])
dataSB <- dataSB %>% mutate_at(vars(matches("Cent")), .funs =function(x) (2*(x-x[1]))/(abs(x)+abs(x[1])))
dataSB <- dataSB %>% mutate_at(vars(matches("Cent")), .funs = function(x) ifelse(is.nan(x), 0, x))
dataSB <- dataSB[-1,] # delete first row of true values
dataSB$Measure <- rep("SB",10)
dataSB$Sample <- 1:10
dataSBLong <- gather(dataSB,key="CentLevel",value="Centrality",starts_with("Cent"), factor_key=TRUE) %>% select(Sample,Network,Measure,Centrality,CentLevel)
```

## Bind all data sets together in long format
```{r DataBind}
hijackPert9 <- bind_rows(dataFTLong, dataFBLong, dataCLLong, dataBWLong, dataSBLong)
hijackPert9$Measure <- ordered(hijackPert9$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.
```{r NPM test}
hijackPert9Lev1 <- subset(hijackPert9, CentLevel == "Cent1", select = c("Sample","Network","Measure","Centrality"))
npmLev1 <- chaAllPairsNashimotoTest(Centrality~Measure,data=hijackPert9Lev1,p.adjust.methods="holm")
hijackPert9Lev2 <- subset(hijackPert9, CentLevel == "Cent2", select = c("Sample","Network","Measure","Centrality"))
npmLev2 <- chaAllPairsNashimotoTest(Centrality~Measure,data=hijackPert9Lev2,p.adjust.methods="holm")
hijackPert9Lev3 <- subset(hijackPert9, CentLevel == "Cent3", select = c("Sample","Network","Measure","Centrality"))
npmLev3 <- chaAllPairsNashimotoTest(Centrality~Measure,data=hijackPert9Lev3,p.adjust.methods="holm")
hijackPert9Lev4 <- subset(hijackPert9, CentLevel == "Cent4", select = c("Sample","Network","Measure","Centrality"))
npmLev4 <- chaAllPairsNashimotoTest(Centrality~Measure,data=hijackPert9Lev4,p.adjust.methods="holm")
hijackPert9Lev5 <- subset(hijackPert9, CentLevel == "Cent5", select = c("Sample","Network","Measure","Centrality"))
npmLev5 <- chaAllPairsNashimotoTest(Centrality~Measure,data=hijackPert9Lev5,p.adjust.methods="holm")
```

# NW Test Results for Level 1
```{r Lev1Results}
npmLev1
```

# NW Test Results for Level 2
```{r Lev2Results}
npmLev2
```

# NW Test Results for Level 3
```{r Lev3Results}
npmLev3
```

# NW Test Results for Level 4
```{r Lev4Results}
npmLev4
```

# NW Test Results for Level 5
```{r Lev5Results}
npmLev5
```


# Line plots to show stability of measures.
Plotting only 5 levels due to overcrowding of plot.
```{r Graphics}
mypal <- c("#CC79A7","#E69F00","#0072B2", "#009E73", "#56B4E9", "#D55E00")
hijackPert9All <- bind_rows(hijackPert9Lev1,hijackPert9Lev2,hijackPert9Lev3,hijackPert9Lev4,hijackPert9Lev5, .id="CentLevel")
ggplot(hijackPert9All, aes(x=Sample,y=Centrality,colour=Measure)) +
  theme_bw() + geom_line(aes(linetype=Measure), size=1.1) + 
  scale_colour_manual(values=mypal) +
  facet_grid(CentLevel~.)
pdf(file="hijackPert9Lines.pdf")
ggplot(hijackPert9All, aes(x=Sample,y=Centrality,colour=Measure)) +
  theme_bw() + geom_line(aes(linetype=Measure), size=1.1) + 
  scale_colour_manual(values=mypal) +
  facet_grid(CentLevel~.)
dev.off()
```

```{r SummaryStats}
sumDataPert9 <- hijackPert9All %>% group_by(CentLevel,Measure) %>% summarise(Mean=mean(Centrality),SD=sd(Centrality),Median=median(Centrality),MAD=mad(Centrality))
ggplot(sumDataPert9,aes(x=CentLevel,y=Mean,colour=Measure)) +
  theme_bw() + geom_point(aes(shape=Measure), size=3) + 
  scale_colour_manual(values=mypal) 
pdf(file="hijackPert9Summary.pdf")
ggplot(sumDataPert9,aes(x=CentLevel,y=Mean,colour=Measure)) +
  theme_bw() + geom_point(aes(shape=Measure), size=3) + 
  scale_colour_manual(values=mypal) 
dev.off()

xtable(sumDataPert9)
```
