---
output:
  html_document: default
  pdf_document: default
---

```{r setup, include=FALSE, echo=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(tidyverse); library(ahp); library(reshape2); library(viridis); library(RColorBrewer)

gg_colors <- function(n) {
		  hues = seq(15, 375, length=n+1)
		  hcl(h=hues, l=65, c=100)[1:n]
		} #color palette

#note set appropriate working directory
#setwd("Analysis")

#############
#############Load and clean data
#############
whales<-read.csv("../Data/Entanglement_data1982-2016", as.is=T)
      x<-as.Date(whales$Observation.Date[1:424], format="%Y/%m/%d")
      x<-c(x, as.Date(whales$Observation.Date[425:nrow(whales)], format="%m/%d/%y"))
      whales$Observation.Date<-x
      whales$Year<-as.numeric(format(whales$Observation.Date, format="%Y"))
      whales$Month<-as.numeric(format(whales$Observation.Date, format="%m"))
      
      whales$Common.Name[grep("humpback", whales$Common.Name, ignore.case=T)]<-"Humpback"
      whales$Common.Name[grep("gray", whales$Common.Name, ignore.case=T)]<-"Gray"	
      whales$Common.Name[grep("blue", whales$Common.Name, ignore.case=T)]<-"Blue"	
      whales$Common.Name[grep("Unidentified", whales$Common.Name, ignore.case=T)]<-"Unidentified"	
      whales$Common.Name[grep("Killer", whales$Common.Name, ignore.case=T)]<-"Killer"
      whales$Common.Name<-gsub(" Whale", "", whales$Common.Name)
      whales$Entanglement.Fishery[grep("dungeness", whales$Entanglement.Fishery, ignore.case=T)]<-"Crab"
      whales$Entanglement.Fishery[grep("DGN", whales$Entanglement.Fishery, ignore.case=T)]<-"Crab"
      whales$Entanglement.Fishery[grep("crab", whales$Entanglement.Fishery, ignore.case=T)]<-"Crab"
      whales$Entanglement.Fishery[grep("sablefish", whales$Entanglement.Fishery, ignore.case=T)]<-"Sablefish"
      whales$Entanglement.Fishery[grep("net", whales$Entanglement.Fishery, ignore.case=T)]<-"Net"
      whales$Entanglement.Fishery[grep("prawn", whales$Entanglement.Fishery, ignore.case=T)]<-"SpotPrawn"
      whales$Entanglement.Fishery[grep("pot", whales$Entanglement.Fishery, ignore.case=T)]<-"Pot"
      whales$Entanglement.Fishery[grep("trap", whales$Entanglement.Fishery, ignore.case=T)]<-"Pot"
      whales$Entanglement.Fishery[grep("unk", whales$Entanglement.Fishery, ignore.case=T)]<-"Unknown"
      whales$Entanglement.Fishery[grep("OTH", whales$Entanglement.Fishery, ignore.case=T)]<-"Unknown"
      whales$Entanglement.Fishery[whales$Entanglement.Fishery==""]<-"Unknown"
      whales$Entanglement.Fishery[whales$Entanglement.Fishery=="HK/LN"]<-"Hook/Line"
      table(whales$Entanglement.Fishery)
      table(whales$Common.Name)
      
      ##Recoding for plotting
      whales$Common.Name[!whales$Common.Name%in%c("Gray", "Humpback", "Unidentified")]<-"Other"  #recode
      whales$Common.Name<-factor(whales$Common.Name, levels=c("Humpback", "Gray", "Other", "Unidentified"))
      whales$Entanglement.Fishery[!whales$Entanglement.Fishery%in%c("Crab", "Net", "Unknown")]<-"Other"  #recode

      
###Load PDO and ENSO data
        enso<-read.table("../Data/MEIindex_fromESRL_NOAA.txt", sep="\t", header=T)
        pdo<-read.table("../Data/PDOindex_http-::www.data.jma.go.jp:gmd:kaiyou:data:db:climate:pdo:pdo.txt", header=F, sep="\t")[,2]
        	pdo<-matrix(pdo, nrow=37, ncol=12, byrow=T) #reshape
          pdo<-as.data.frame(pdo); row.names(pdo) <- seq(1980, 2016, 1); colnames(pdo) <- month.name
	
```


Wilcoxon test to measure degree of increase in Humpback entanglements, pre- vs -post-2012

```{r}
wilcox.test(
table(format(whales$Observation.Date[whales$Common.Name=="Humpback"], "%Y"))[1:19],
table(format(whales$Observation.Date[whales$Common.Name=="Humpback"], "%Y"))[20:22]
)

```

Annual reported Humpback entanglements have increased significantly in recent years (1984-2012 vs. 2014-2106; Wilcoxon test; p = .007)


Figure 1: Trend over time
```{r, echo=FALSE}
#OVERALL
pdf("../Figures/Figure1_Species_Year.pdf", width = 8, height = 4)
ggplot(whales,aes(x= Year, fill=Common.Name))+
  geom_histogram(binwidth=2, position = "dodge")+
  scale_fill_manual(values = RColorBrewer::brewer.pal(9, "Blues")[c(9, 7, 5, 3)], name = "Species")+
  theme_bw()+
  #theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), 
  #      panel.background = element_blank(), axis.line = element_line(colour = "black"))+
  theme(axis.text.x = element_text(angle = 90, hjust = 1))+
  theme(strip.text.x = element_blank())+
  labs(x = "Year", y = "Reported Entanglements (Count)")
dev.off()
```


Figure 2: Trend by time and by fishery
```{r, echo=FALSE}
pdf("../Figures/Figure2_Species_Year_Fishery.pdf", width = 5, height = 4)
ggplot(whales[whales$Common.Name%in%c("Humpback","Gray"),],aes(x= Year, fill=Common.Name))+
  geom_histogram(binwidth=1, position = "dodge")+
  scale_fill_manual(values = RColorBrewer::brewer.pal(9, "Blues")[c(9, 7, 5, 3)], name = "Species")+
  facet_grid(Entanglement.Fishery~Common.Name)+
  theme_bw()+
  #theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), 
  #      panel.background = element_blank(), axis.line = element_line(colour = "black"))+
  theme(axis.text.x = element_text(angle = 90, hjust = 1))+
  theme(strip.text.x = element_blank())+
  labs(x = "Year", y = "Reported Entanglements (Count)")
dev.off()

```
\vspace{2em}


Supplemental figure:
The increase over previous years is most striking in May, June, August, and September. 
```{r, echo=FALSE}
MONTHS <- month.name; names(MONTHS)<-seq(1,12,1)

pdf("../Supplemental_Info/SupplFig1.pdf", width = 8, height = 6)
   ggplot(whales[whales$Common.Name%in%c("Humpback", "Gray"),],aes(x= Year, fill=Common.Name))+
     geom_histogram(binwidth=2, position = "dodge")+
     facet_wrap(~Month, labeller = as_labeller(MONTHS))+
     theme_bw()+
     theme(axis.text.x = element_text(angle = 90, hjust = 1))+
     scale_fill_manual(values = RColorBrewer::brewer.pal(9, "Blues")[c(9, 5)], name = "Species")+
     labs(x = "Year", y = "Reported Entanglements (Count)")
dev.off()


```
\vspace{2em}

No real relationship between entanglements and PDO or ENSO in any month of the year.  3 outlier years (the recent ones) might cause apparent correlation, but this doesn't hold up when looking at the rest of the data.

```{r}
HB<-as.data.frame(table(whales[whales$Common.Name=="Humpback","Year"]))
  names(HB)<-c("YEAR","Count")

  
###ENSO Regresssion
  HB_enso<-merge(HB, enso)
    summary(lm(log(HB_enso[,2])~HB_enso[,3]))$coefficients[8]
  ensoRegression <- lapply(HB_enso[,3:14], FUN = function(x) summary(lm(log(HB_enso[,2])~x))$coefficients[8])  #get p-value for linear model of (log) whale entanglements (dependent variable) vs. ENSO index for each month (independent variable).  Result is 12 values
  
  summary(unlist(ensoRegression))  #no significant correlation
  
  #Plot this, if desired
  HB_enso<-melt(HB_enso, id.vars = c("YEAR", "Count"))
  HB_enso$YEAR<-as.numeric(as.character(HB_enso$YEAR))
  ggplot(HB_enso, aes(x=Count, y=value))+geom_point()+
    geom_smooth(method = "lm")+
    scale_x_continuous(trans='log10')+
    facet_wrap(~variable)

###PDO Regresssion  
  
  pdo<-data.frame(YEAR = row.names(pdo), pdo)
  pdo<-merge(HB, pdo)

  pdoRegression <- lapply(pdo[,3:14], FUN = function(x) summary(lm(log(pdo[,2])~x))$coefficients[8])  #get p-value for linear model of (log) whale entanglements (dependent variable) vs. pdo index for each month (independent variable).  Result is 12 values
  summary(unlist(pdoRegression)) #no significant correlation
  
  #Plot this, if desired
    pdo<-melt(pdo, id.vars = c("YEAR", "Count"))
    pdo$YEAR<-as.numeric(as.character(pdo$YEAR))

    ggplot(pdo, aes(x=Count, y=value))+geom_point()+
    geom_smooth(method = "lm")+
    scale_x_continuous(trans='log10')+
    facet_wrap(~variable)

  
```



###MCDA analysis using package "ahp"

```{r}
##Evaluate policy alternatives under many different weighting schemes

templateFile<-readLines("MCDA/WhaleEntanglement_20190329.ahp") #read in template file for MCDA analysis using package "ahp"
  ChangeWeights <- function(y, x) gsub("1", y, templateFile[x])  #create function to change weights in that file


set.seed(49)
resultsList<-list(AssignedWeights = list(), ProportionalWeights = list(), RankedAlternatives = list()) #create empty list to store results

for (k in 1:100){  
    tempFile<-templateFile  
    weights<-round(rlnorm(3),1) #sample from a lognormal distribution, so that weight for a given criterion cannot be zero
        #c(1,1,5) #a vector of 3 values representing the relative weights for cost-vs-effectiveness, cost-vs-response, and effectiveness-vs-response, respectively. Following the "ahp" syntax,     
        # preferences are defined pairwise, where:
        # 1 means: A is equal to B
        # 9 means: A is highly preferrable to B
        # 1/9 means: B is highly preferrable to A
    
        lines<-99:101  #adjust weights as appropriate  NOTE: these line numbers change when you add/subtract alternatives, etc.  TODO: use grep
          for (i in 1:3){
            tempFile[lines[i]] <- ChangeWeights(weights[i], lines[i])
        }
    
        #Conduct MCDA analysis using package "ahp"
        tempAnalysis<-LoadString(tempFile)
        Calculate(tempAnalysis)
        tempResults<-as.data.frame(AnalyzeTable(tempAnalysis, sort = "orig")) #keep order of alternatives as in the input data
        
        #store in resultsList
        resultsList[["AssignedWeights"]][[k]]<-weights
        resultsList$ProportionalWeights[[k]]<-tempResults[2:4,2]  #in order of Cost, Effectiveness, Response
        resultsList$RankedAlternatives[[k]]<-order(tempResults[1,3:18], decreasing = T) #preserve the rank order of alternatives under the weighting scenario tested
}

```


```{r, warning=F, echo = F, message=F}
####Plot the results showing ranks of policy alternatives under scenarios of different weightings


weights.df<-as.data.frame(resultsList$ProportionalWeights)
  row.names(weights.df)<-c("Cost", "Effectiveness", "Response")
  weights.df<-weights.df[,order(weights.df[1,], weights.df[2,], decreasing = T)] #order first by cost, then by effectiveness
  
  
rank.df<-as.data.frame(resultsList$RankedAlternatives)
  row.names(rank.df)<-names(tempResults)[3:18]
  rank.df<-rank.df[,order(weights.df[1,], weights.df[2,], decreasing = T)] #reorder to match the order of the weighting scenarios
  colnames(rank.df)<-paste("WeightScenario", seq(1,100,1), sep = "_")
  
  labelAlternatives <- sapply(1:nrow(rank.df), function(x) strsplit(row.names(rank.df), " ")[[x]][2])
```  

```{r FrequencyRankTable}

#how often does each policy option rank #1, #2, etc, in order of preference?
  
freq.rank.table <- as.data.frame(matrix(NA, nrow = 16, ncol = 15))
  for (i in 1:16){
    for(j in 1:16){
      freq.rank.table[i,j] <- sum(rank.df[i,] == j)
    }
  }
  row.names(freq.rank.table) <- labelAlternatives
  colnames(freq.rank.table) <- paste(sprintf("%02d", 1:16))  #note leading 0 in numbering allows for easy sorting
  
  
```


Figure 3: MCDA Results and Weighting Schemes

```{r}
  colorScale <- gg_colors(16)
    #c("#65cfd8", "#eaa99e", "#8be9d4", "#d9b2e9", "#b3dba3", "#ebaec8", "#8bc5ac", "#aabbea", "#ede1a3", "#7cccee", "#e8b893", "#b3e9db", "#bebc81", "#e3e2b9", "#c3bc96")  #calculated using http://tools.medialab.sciences-po.fr/iwanthue/ , fancy presets, 15 colors


pdf("../Figures/Figure3_MCDA.pdf", width = 6.3, height = 5)
    nf <- layout(matrix(c(1,1,1,
                          1,1,1,
                          2,2,2), 
                        3, 3, byrow = TRUE), respect = TRUE)
    
    #Main results Plot
              par(mar=c(0,2,0,6), 
                  mgp = c(1,-1,-1), 
                  lab=c(5,16,7)
              )
              
              #vector for line weights 
              lwdvec <- 5/log(rowSums(rank.df))
              
              plot(as.numeric(rank.df[1,]), type="l", ylim=c(16,1), 
                   col= colorScale[1], #viridis(15, alpha = .7)[1], 
                   bty="n", ylab="Rank", xlab="", xaxt="n", yaxt="n", lwd = lwdvec[1])
              axis(2, lwd = 0, lwd.ticks = 0)
              mtext("A", side = 3, line = -2.5, at = -7, cex = 2)
              
              
              for (i in 2:16){
                points(as.numeric(rank.df[i,]), type="l", col=colorScale[i], lwd = lwdvec[i])  #viridis(15, alpha = .7)[i]
              }
              legend(x=100, y=3.5, legend = labelAlternatives , fill = colorScale, bty="n", ncol=1, xpd=NA)
              #par(xpd=TRUE)
              #text(x = 135, y=8.5, labels = "Policy Alternatives", srt = 270)
              
    
    #Plot of weights as the x-axis label
    
              par(mar=c(0,2,0,6))
              barplot(as.matrix(weights.df), xlab="", ylab = "Weight", 
                      xaxt = "n", yaxt = "n", bty="L", 
                      col = colorScale[c(9,15,11)], #viridis(15, alpha = .7)[c(1,5,9)], 
                      border = NA, space = 0)
              mtext("B", side = 3, line = -2.5, at = -7, cex = 2)
              
              legend(x = 100, y=.5, legend = c("Cost", "Effectiveness", "Response"), fill = colorScale[c(9,15,11)], #viridis(15, alpha = .7)[c(1,5,10)], 
                     bty="n", ncol=1, xpd=NA)
dev.off()
```

Figure 4: Frequency Rank Table for Policy Alternatives
  
```{r Plotting}

pdf("../Figures/Figure4_AlternativesRankFrequency.pdf", width = 6.3, height = 5)
freq.rank.table %>% 
    rownames_to_column(var = "PolicyAlternative") %>% 
    gather(key = "Rank", value = "Frequency", -PolicyAlternative) %>% 
    mutate(Rank = as.numeric(Rank)) %>% 
    arrange(desc(Frequency * (1/Rank))) %>% 
    ggplot(aes(y= Rank, #factor(Rank, levels = as.character(unique(Rank))),
               x = factor(PolicyAlternative, levels = unique(PolicyAlternative)), 
               fill = Frequency)) +
      geom_tile() +
      scale_fill_gradient(low = "#EFF3FF", high = "#2171B5") +
      ylab("Rank") + xlab("Policy Alternative") +
      theme_bw() +
      theme(axis.text.x = element_text(angle = 90, hjust = 1),
            panel.grid.major = element_blank(),
            panel.grid.minor = element_blank(),
            panel.border = element_blank())

dev.off()
```
















