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”.

# setwd("./Code/HijackerData")
#install.packages("coin") - use if you have not installed the package.
library(coin) # load the package in R 
## Loading required package: survival
library(tidyverse) # for tidying the data
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ ggplot2 3.3.5     ✓ purrr   0.3.4
## ✓ tibble  3.1.5     ✓ dplyr   1.0.7
## ✓ tidyr   1.1.4     ✓ stringr 1.4.0
## ✓ readr   2.0.2     ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
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 12% perturbation and calculate relative percentage error for each centrality measure.

Create variable names for all data sets

varNames <- c("Network",paste0(c("Cent","Tier"), rep(1:62,each=2)))

Import betweenness data and convert to long format

dataBW <- read.csv(file="Hijack_12percent_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])
## [1] 11
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

dataCL <- read.csv(file="Hijack_12percent_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])
## [1] 11
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

dataFB <- read.csv(file="Hijack_12percent_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])
## [1] 11
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

dataFT <- read.csv(file="Hijack_12percent_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])
## [1] 11
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

dataSB <- read.csv(file="Hijack_12percent_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])
## [1] 92
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

hijackPert12 <- bind_rows(dataFTLong, dataFBLong, dataCLLong, dataBWLong, dataSBLong)
hijackPert12$Measure <- ordered(hijackPert12$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.

hijackPert12Lev1 <- subset(hijackPert12, CentLevel == "Cent1", select = c("Sample","Network","Measure","Centrality"))
npmLev1 <- chaAllPairsNashimotoTest(Centrality~Measure,data=hijackPert12Lev1,p.adjust.methods="holm")
hijackPert12Lev2 <- subset(hijackPert12, CentLevel == "Cent2", select = c("Sample","Network","Measure","Centrality"))
npmLev2 <- chaAllPairsNashimotoTest(Centrality~Measure,data=hijackPert12Lev2,p.adjust.methods="holm")
hijackPert12Lev3 <- subset(hijackPert12, CentLevel == "Cent3", select = c("Sample","Network","Measure","Centrality"))
npmLev3 <- chaAllPairsNashimotoTest(Centrality~Measure,data=hijackPert12Lev3,p.adjust.methods="holm")
hijackPert12Lev4 <- subset(hijackPert12, CentLevel == "Cent4", select = c("Sample","Network","Measure","Centrality"))
npmLev4 <- chaAllPairsNashimotoTest(Centrality~Measure,data=hijackPert12Lev4,p.adjust.methods="holm")
hijackPert12Lev5 <- subset(hijackPert12, CentLevel == "Cent5", select = c("Sample","Network","Measure","Centrality"))
npmLev5 <- chaAllPairsNashimotoTest(Centrality~Measure,data=hijackPert12Lev5,p.adjust.methods="holm")

NW Test Results for Level 1

npmLev1
## 
##  Pairwise comparisons using Pairwise comparisons using Nashimoto-Wright's NPT'-Test
## data: Centrality by Measure
##    FT      FB      CL      BW     
## FB 1.00000 -       -       -      
## CL 1.00000 1.00000 -       -      
## BW 0.00063 0.00063 0.00063 -      
## SB 0.00063 0.00063 0.00063 1.00000
## 
## P value adjustment method: holm
## alternative hypothesis: greater

NW Test Results for Level 2

npmLev2
## 
##  Pairwise comparisons using Pairwise comparisons using Nashimoto-Wright's NPT'-Test
## data: Centrality by Measure
##    FT      FB      CL      BW   
## FB 0.779   -       -       -    
## CL 1.000   1.000   -       -    
## BW 0.082   0.082   0.082   -    
## SB 1.3e-05 1.3e-05 1.3e-05 0.045
## 
## P value adjustment method: holm
## alternative hypothesis: greater

NW Test Results for Level 3

npmLev3
## 
##  Pairwise comparisons using Pairwise comparisons using Nashimoto-Wright's NPT'-Test
## data: Centrality by Measure
##    FT      FB      CL      BW     
## FB 1.00000 -       -       -      
## CL 1.00000 1.00000 -       -      
## BW 0.00091 0.00091 0.00091 -      
## SB 0.02214 0.02214 0.02214 1.00000
## 
## P value adjustment method: holm
## alternative hypothesis: greater

NW Test Results for Level 4

npmLev4
## 
##  Pairwise comparisons using Pairwise comparisons using Nashimoto-Wright's NPT'-Test
## data: Centrality by Measure
##    FT   FB   CL   BW  
## FB 1.00 -    -    -   
## CL 1.00 1.00 -    -   
## BW 0.48 0.48 0.48 -   
## SB 0.11 0.11 0.11 0.83
## 
## P value adjustment method: holm
## alternative hypothesis: greater

NW Test Results for Level 5

npmLev5
## 
##  Pairwise comparisons using Pairwise comparisons using Nashimoto-Wright's NPT'-Test
## data: Centrality by Measure
##    FT    FB    CL    BW   
## FB 1.000 -     -     -    
## CL 1.000 1.000 -     -    
## BW 0.808 0.808 1.000 -    
## SB 0.035 0.035 0.118 0.387
## 
## P value adjustment method: holm
## alternative hypothesis: greater

Line plots to show stability of measures.

Plotting only 5 levels due to overcrowding of plot.

mypal <- c("#CC79A7","#E69F00","#0072B2", "#009E73", "#56B4E9", "#D55E00")
hijackPert12All <- bind_rows(hijackPert12Lev1,hijackPert12Lev2,hijackPert12Lev3,hijackPert12Lev4,hijackPert12Lev5, .id="CentLevel")
ggplot(hijackPert12All, 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="hijackPert12Lines.pdf")
ggplot(hijackPert12All, 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()
## quartz_off_screen 
##                 2
sumDataPert12 <- hijackPert12All %>% group_by(CentLevel,Measure) %>% summarise(Mean=mean(Centrality),SD=sd(Centrality),Median=median(Centrality),MAD=mad(Centrality))
## `summarise()` has grouped output by 'CentLevel'. You can override using the `.groups` argument.
ggplot(sumDataPert12,aes(x=CentLevel,y=Mean,colour=Measure)) +
  theme_bw() + geom_point(aes(shape=Measure), size=3) + 
  scale_colour_manual(values=mypal) 
## Warning: Using shapes for an ordinal variable is not advised

pdf(file="hijackPert12Summary.pdf")
ggplot(sumDataPert12,aes(x=CentLevel,y=Mean,colour=Measure)) +
  theme_bw() + geom_point(aes(shape=Measure), size=3) + 
  scale_colour_manual(values=mypal) 
## Warning: Using shapes for an ordinal variable is not advised
dev.off()
## quartz_off_screen 
##                 2
xtable(sumDataPert12)
## % latex table generated in R 4.1.2 by xtable 1.8-4 package
## % Tue Jul 12 12:04:03 2022
## \begin{table}[ht]
## \centering
## \begin{tabular}{rllrrrr}
##   \hline
##  & CentLevel & Measure & Mean & SD & Median & MAD \\ 
##   \hline
## 1 & 1 & FT & 0.00 & 0.00 & 0.00 & 0.00 \\ 
##   2 & 1 & FB & 0.00 & 0.00 & 0.00 & 0.00 \\ 
##   3 & 1 & CL & -0.27 & 0.61 & -0.05 & 0.03 \\ 
##   4 & 1 & BW & 0.00 & 0.00 & 0.00 & 0.00 \\ 
##   5 & 1 & SB & 0.00 & 0.00 & 0.00 & 0.00 \\ 
##   6 & 2 & FT & 0.09 & 0.14 & 0.09 & 0.14 \\ 
##   7 & 2 & FB & 0.16 & 0.23 & 0.21 & 0.14 \\ 
##   8 & 2 & CL & -0.28 & 0.61 & -0.06 & 0.06 \\ 
##   9 & 2 & BW & 0.14 & 0.29 & 0.17 & 0.36 \\ 
##   10 & 2 & SB & 0.68 & 0.53 & 0.58 & 0.44 \\ 
##   11 & 3 & FT & 0.60 & 0.97 & 0.00 & 0.00 \\ 
##   12 & 3 & FB & -0.08 & 0.42 & 0.02 & 0.08 \\ 
##   13 & 3 & CL & -0.29 & 0.61 & -0.06 & 0.06 \\ 
##   14 & 3 & BW & 0.60 & 0.97 & 0.00 & 0.00 \\ 
##   15 & 3 & SB & 0.00 & 0.00 & 0.00 & 0.00 \\ 
##   16 & 4 & FT & 0.03 & 0.18 & 0.03 & 0.19 \\ 
##   17 & 4 & FB & 0.09 & 0.31 & 0.06 & 0.36 \\ 
##   18 & 4 & CL & -0.26 & 0.61 & -0.04 & 0.02 \\ 
##   19 & 4 & BW & 0.22 & 0.66 & 0.01 & 0.48 \\ 
##   20 & 4 & SB & 0.31 & 0.78 & 0.14 & 0.52 \\ 
##   21 & 5 & FT & -0.49 & 1.06 & -0.09 & 0.69 \\ 
##   22 & 5 & FB & -0.94 & 0.87 & -1.24 & 0.85 \\ 
##   23 & 5 & CL & -0.25 & 0.62 & -0.04 & 0.04 \\ 
##   24 & 5 & BW & -0.51 & 1.38 & -0.14 & 1.99 \\ 
##   25 & 5 & SB & 0.00 & 0.00 & 0.00 & 0.00 \\ 
##    \hline
## \end{tabular}
## \end{table}