## Entering flow centrality simulations
## Files are in "/Users/Monnie/Dropbox/2018 Fall Leave/FlowThroughCentrality/Simulations"
## Data are in https://smu.box.com/s/zjktbic32bnefnb9iv4wk3pkceubeveu
library(tidyverse)

setwd("../Nodes50Deg6Pert6")
# Enter all data files
varNames <- c("Network",paste0(c("Cent","Tier"), rep(1:50,each=2)))
dataBW <- read.csv(file="Nodes50_AvgDeg06_UW_01_Pert06_Betweennesss.csv",skip=5,nrows=100,header=FALSE)
colnames(dataBW) <- varNames
dataBW$Measure <- rep("BW",100)
dataCL <- read.csv(file="Nodes50_AvgDeg06_UW_01_Pert06_Closeness.csv",skip=5,nrows=100,header=FALSE)
colnames(dataCL) <- varNames
dataCL$Measure <- rep("CL",100)
dataFB <- read.csv(file="Nodes50_AvgDeg06_UW_01_Pert06_FlowBetween.csv",skip=5,nrows=100,header=FALSE)
colnames(dataFB) <- varNames
dataFB$Measure <- rep("FB",100)
dataFT <- read.csv(file="Nodes50_AvgDeg06_UW_01_Pert06_Flowthrough.csv",skip=5,nrows=100,header=FALSE)
colnames(dataFT) <- varNames
dataFT$Measure <- rep("FT",100)
dataSB <- read.csv(file="Nodes50_AvgDeg06_UW_01_Pert06_StableBetween.csv",skip=5,nrows=100,header=FALSE)
colnames(dataSB) <- varNames
dataSB$Measure <- rep("SB",100)

## Create tibbles with first three columns of each file
dataBWfirst3 <- dplyr::select(dataBW,c(Network:Tier3,Measure))
dataCLfirst3 <- dplyr::select(dataCL,c(Network:Tier3,Measure))
dataFBfirst3 <- dplyr::select(dataFB,c(Network:Tier3,Measure))
dataFTfirst3 <- dplyr::select(dataFT,c(Network:Tier3,Measure))
dataSBfirst3 <- dplyr::select(dataSB,c(Network:Tier3,Measure))
dataAllFirst3 <- bind_rows(dataBWfirst3,dataCLfirst3,dataFBfirst3,dataFTfirst3, dataSBfirst3)
dataAllFirst3$Measure <- as.factor(dataAllFirst3$Measure)

## Data in long format
tier3n50deg6pert6a <- gather(dataAllFirst3,key=CentLevel,value=Centrality,Cent1,Cent2,Cent3,factor_key=TRUE) 
tier3n50deg6pert6b <- gather(dataAllFirst3,key=Level,value=Tier, Tier1,Tier2,Tier3, factor_key=TRUE)
tier3n50deg6pert6 <-merge(tier3n50deg6pert6a,tier3n50deg6pert6b,by=c("Network","Measure"))
tier3n50deg6pert6 <- dplyr::select(tier3n50deg6pert6,-c(Tier1:Tier3,Cent1:Cent3))
tier3n50deg6pert6 <- separate(tier3n50deg6pert6,Network,into=c("Trash","Sample"),sep=" ")
tier3n50deg6pert6$Sample <- as.integer(tier3n50deg6pert6$Sample)
tier3n50deg6pert6 <- tier3n50deg6pert6[,-1]

pdf(file="tier3n50deg6pert6Lines.pdf")
ggplot(tier3n50deg6pert6,aes(x=Sample, y=Centrality, colour=CentLevel)) + geom_line() + facet_grid(Measure~.,scales="free")
dev.off()

xtable(tier3n50deg6pert6 %>% group_by(CentLevel,Measure) %>% summarise(mean=mean(Centrality),sd=sd(Centrality),cv=100*(sd(Centrality)/mean(Centrality))))

## More plots
http://www.sthda.com/english/wiki/ggplot2-error-bars-quick-start-guide-r-software-and-data-visualization

## Some histograms. Data are not Normal for the most part.
par(mfrow=c(5,3))
hist(dataAllFirst3$BWTier1)
hist(dataAllFirst3$BWTier2)
hist(dataAllFirst3$BWTier3)
hist(dataAllFirst3$FBTier1)
hist(dataAllFirst3$FBTier2)
hist(dataAllFirst3$FBTier3)
hist(dataAllFirst3$FTTier1)
hist(dataAllFirst3$FTTier2)
hist(dataAllFirst3$FTTier3)
hist(dataAllFirst3$SBTier1)
hist(dataAllFirst3$SBTier2)
hist(dataAllFirst3$SBTier3)
hist(dataAllFirst3$CLTier1)
hist(dataAllFirst3$CLTier2)
hist(dataAllFirst3$CLTier3)

## Create tibble in long format

## Run tests for variance on dependent variables

## References
## http://psychologicalstatistics.blogspot.com/2006/05/what-is-all-this-stuff-about.html
### https://www.rdocumentation.org/packages/vegan/versions/2.4-2/topics/mrpp
## http://cc.oulu.fi/~jarioksa/softhelp/vegan/html/mrpp.html
