#########################################
#                                       #
#   Mean trophic level function         #
#                                       #
#        created:  20/12/2016           #
#          LANDINGS data_FR             #
#  All species & no Pelagic species     #
#                                       #
#                                       #
#########################################

#! FANCY AXIS FUNCTION
"fancyAxis" <- function(side, at, skip = 2, ...)
{
  axis(side, at = at, labels = FALSE)
  axis(side, at = at[ ( (1:length(at) %% skip ) == 1 )],
       labels = TRUE, tick = FALSE, ...)
}


TL<-read.table("C:/[...]/TL_landings&EVHOE_FR.csv",sep=";",dec=",",header=T)
biomass2<-read.table("C:/[...]/ICES_landings1997_2014_AllSpecies_FR.csv",sep=";",dec=",",header=T)
biomass3<-read.table("C:/[...]/ICES_landings1997_2014_noPelagic_FR.csv",sep=";",dec=",",header=T)
View(TL)
View(biomass2)
View(biomass3)


############################################################################################
###################  MTL_2.0 Landings   ####################################
############################################################################################
cutoff2<-2.0
rep<-500

biomass2$year<-as.factor(biomass2$year)
rownames(TL)<-TL$species
year<-unique(biomass2$year)

MTL2land<-matrix(nrow=18,ncol=rep,NA)
MTL2land<-data.frame(MTL2land)
row.names(MTL2land)<-year
for (k in 1:rep)
  for (i in year)
  { print(paste("rep: ",k,"    year: ",i,sep=""))
    selection <- biomass2[which(biomass2$year==i),]
    row.names(selection) <- selection$species
    TLselect <- TL[TL$TL >= cutoff2, ]
    selection <- selection[selection$species %in% TLselect$species, ]
    
    biomass_year<-sum(selection$biomass)
    for (j in 1:dim(selection)[1])
    {
      SP<-selection[j,"species"]
      TL_sel<-TL[paste(SP),"TL"]
      SE_sel<-TL[paste(SP),"SE"]
      TL_sample<-sample(seq(TL_sel-SE_sel,TL_sel+SE_sel,0.0001),1)
      selection[j,"biomass"]<-(selection[j,"biomass"]*TL_sample)
    }
    MTL2land[i,k]<-(sum(selection[,"biomass"])/biomass_year)
  }

############################################################################################
###################  MTL_3.25  Landings ####################################
############################################################################################
cutoff3<-3.25
rep<-500

biomass2$year<-as.factor(biomass2$year)
rownames(TL)<-TL$species
year<-unique(biomass2$year)

MTL3land<-matrix(nrow=18,ncol=rep,NA)
MTL3land<-data.frame(MTL3land)
row.names(MTL3land)<-year
for (k in 1:rep)
  for (i in year)
  { print(paste("rep: ",k,"    year: ",i,sep=""))
    selection <- biomass2[which(biomass2$year==i),]
    row.names(selection) <- selection$species
    TLselect <- TL[TL$TL >= cutoff3, ]
    selection <- selection[selection$species %in% TLselect$species, ]
    
    biomass_year<-sum(selection$biomass)
    for (j in 1:dim(selection)[1])
    {
      SP<-selection[j,"species"]
      TL_sel<-TL[paste(SP),"TL"]
      SE_sel<-TL[paste(SP),"SE"]
      TL_sample<-sample(seq(TL_sel-SE_sel,TL_sel+SE_sel,0.0001),1)
      selection[j,"biomass"]<-(selection[j,"biomass"]*TL_sample)
    }
    MTL3land[i,k]<-(sum(selection[,"biomass"])/biomass_year)
  }

############################################################################################
###################  MTL_4   Landings ####################################
############################################################################################
cutoff4<-4.0
rep<-500

biomass2$year<-as.factor(biomass2$year)
rownames(TL)<-TL$species
year<-unique(biomass2$year)

MTL4land<-matrix(nrow=18,ncol=rep,NA)
MTL4land<-data.frame(MTL4land)
row.names(MTL4land)<-year
for (k in 1:rep)
  for (i in year)
  { print(paste("rep: ",k,"    year: ",i,sep=""))
    selection <- biomass2[which(biomass2$year==i),]
    row.names(selection) <- selection$species
    TLselect <- TL[TL$TL >= cutoff4, ]
    selection <- selection[selection$species %in% TLselect$species, ]
    
    biomass_year<-sum(selection$biomass)
    for (j in 1:dim(selection)[1])
    {
      SP<-selection[j,"species"]
      TL_sel<-TL[paste(SP),"TL"]
      SE_sel<-TL[paste(SP),"SE"]
      TL_sample<-sample(seq(TL_sel-SE_sel,TL_sel+SE_sel,0.0001),1)
      selection[j,"biomass"]<-(selection[j,"biomass"]*TL_sample)
    }
    MTL4land[i,k]<-(sum(selection[,"biomass"])/biomass_year)
  }

############################################################################################
###################  MTL_2.0 Landings no Pelagic   ####################################
############################################################################################

biomass3$year<-as.factor(biomass3$year)
rownames(TL)<-TL$species
year<-unique(biomass3$year)

MTL2landnoP<-matrix(nrow=18,ncol=rep,NA)
MTL2landnoP<-data.frame(MTL2landnoP)
row.names(MTL2landnoP)<-year
for (k in 1:rep)
  for (i in year)
  { print(paste("rep: ",k,"    year: ",i,sep=""))
    selection <- biomass3[which(biomass3$year==i),]
    row.names(selection) <- selection$species
    TLselect <- TL[TL$TL >= cutoff2, ]
    selection <- selection[selection$species %in% TLselect$species, ]
    
    biomass_year<-sum(selection$biomass)
    for (j in 1:dim(selection)[1])
    {
      SP<-selection[j,"species"]
      TL_sel<-TL[paste(SP),"TL"]
      SE_sel<-TL[paste(SP),"SE"]
      TL_sample<-sample(seq(TL_sel-SE_sel,TL_sel+SE_sel,0.0001),1)
      selection[j,"biomass"]<-(selection[j,"biomass"]*TL_sample)
    }
    MTL2landnoP[i,k]<-(sum(selection[,"biomass"])/biomass_year)
  }

############################################################################################
###################  MTL_3.25  Landings no Pelagic ####################################
############################################################################################

biomass3$year<-as.factor(biomass3$year)
rownames(TL)<-TL$species
year<-unique(biomass3$year)

MTL3landnoP<-matrix(nrow=18,ncol=rep,NA)
MTL3landnoP<-data.frame(MTL3landnoP)
row.names(MTL3landnoP)<-year
for (k in 1:rep)
  for (i in year)
  { print(paste("rep: ",k,"    year: ",i,sep=""))
    selection <- biomass3[which(biomass3$year==i),]
    row.names(selection) <- selection$species
    TLselect <- TL[TL$TL >= cutoff3, ]
    selection <- selection[selection$species %in% TLselect$species, ]
    
    biomass_year<-sum(selection$biomass)
    for (j in 1:dim(selection)[1])
    {
      SP<-selection[j,"species"]
      TL_sel<-TL[paste(SP),"TL"]
      SE_sel<-TL[paste(SP),"SE"]
      TL_sample<-sample(seq(TL_sel-SE_sel,TL_sel+SE_sel,0.0001),1)
      selection[j,"biomass"]<-(selection[j,"biomass"]*TL_sample)
    }
    MTL3landnoP[i,k]<-(sum(selection[,"biomass"])/biomass_year)
  }

############################################################################################
###################  MTL_4   Landings no Pelagic ####################################
############################################################################################

biomass3$year<-as.factor(biomass3$year)
rownames(TL)<-TL$species
year<-unique(biomass3$year)

MTL4landnoP<-matrix(nrow=18,ncol=rep,NA)
MTL4landnoP<-data.frame(MTL4landnoP)
row.names(MTL4landnoP)<-year
for (k in 1:rep)
  for (i in year)
  { print(paste("rep: ",k,"    year: ",i,sep=""))
    selection <- biomass3[which(biomass3$year==i),]
    row.names(selection) <- selection$species
    TLselect <- TL[TL$TL >= cutoff4, ]
    selection <- selection[selection$species %in% TLselect$species, ]
    
    biomass_year<-sum(selection$biomass)
    for (j in 1:dim(selection)[1])
    {
      SP<-selection[j,"species"]
      TL_sel<-TL[paste(SP),"TL"]
      SE_sel<-TL[paste(SP),"SE"]
      TL_sample<-sample(seq(TL_sel-SE_sel,TL_sel+SE_sel,0.0001),1)
      selection[j,"biomass"]<-(selection[j,"biomass"]*TL_sample)
    }
    MTL4landnoP[i,k]<-(sum(selection[,"biomass"])/biomass_year)
  }


###################################################
########### PLOT (after calculations) #############
###################################################

####################
# PLOT MTL_2.0    #
####################

years2land <- as.numeric(rownames(MTL2land))
MTL2land[is.na(MTL2land)] <- 4

par(mar=c(6,7,2,2), mgp=c(4.5, 1, 0), las=0)

plot(NULL, xlim = c(min(years2land), max(years2land)), ylim=c(3.5,4.7),
     ylab="MTL_2.0",xlab="Year",cex.lab = 3,
     axes=FALSE, frame.plot = TRUE, bty = "l")

stdevs2land <- apply(MTL2land, 1, sd)
lcl2land <- rowMeans(MTL2land) - stdevs2land
ucl2land <- rowMeans(MTL2land) + stdevs2land
x.polygon2land <- c(years2land, sort(years2land, decreasing = TRUE))
y.polygon2land <- c(ucl2land, rev(lcl2land))

polygon(x.polygon2land, y.polygon2land, col = "#F2E399", border = NA)
lines(years2land, rowMeans(MTL2land), col = "#006349", lwd = 3)

fancyAxis(1, at = seq(1997, 2015, 2), las = 1, cex.axis = 2.5)
fancyAxis(2, at = seq(3.5, 4.7, 0.2), las = 1, cex.axis = 2.5)

rect(1998, 3.5, 2000, 4.7, col = "white", border = NA)


####################
# PLOT MTL_3.25    #
####################

years3land <- as.numeric(rownames(MTL3land))
MTL3land[is.na(MTL3land)] <- 4

par(mar=c(6,7,2,2), mgp=c(4.5, 1, 0), las=0)

plot(NULL, xlim = c(min(years3land), max(years3land)), ylim=c(3.5,4.7),
     ylab="MTL_3.25",xlab="Year",cex.lab = 3,
     axes=FALSE, frame.plot = TRUE, bty = "l")

stdevs3land <- apply(MTL3land, 1, sd)
lcl3land <- rowMeans(MTL3land) - stdevs3land
ucl3land <- rowMeans(MTL3land) + stdevs3land
x.polygon3land <- c(years3land, sort(years3land, decreasing = TRUE))
y.polygon3land <- c(ucl3land, rev(lcl3land))

polygon(x.polygon3land, y.polygon3land, col = "#F2E399", border = NA)
lines(years3land, rowMeans(MTL3land), col = "#006349", lwd = 3)

fancyAxis(1, at = seq(1997, 2015, 2), las = 1, cex.axis = 2.5)
fancyAxis(2, at = seq(3.5, 4.7, 0.2), las = 1, cex.axis = 2.5)

rect(1998, 3.5, 2000, 4.7, col = "white", border = NA)


####################
# PLOT MTL_4.0    #
####################

years4land <- as.numeric(rownames(MTL4land))
MTL4land[is.na(MTL4land)] <- 4.2

par(mar=c(6,7,2,2), mgp=c(4.5, 1, 0), las=0)

plot(NULL, xlim = c(min(years4land), max(years4land)), ylim=c(3.5,4.7),
     ylab="MTL_4.0",xlab="Year",cex.lab = 3,
     axes=FALSE, frame.plot = TRUE, bty = "l")

stdevs4land <- apply(MTL4land, 1, sd)
lcl4land <- rowMeans(MTL4land) - stdevs4land
ucl4land <- rowMeans(MTL4land) + stdevs4land
x.polygon4land <- c(years4land, sort(years4land, decreasing = TRUE))
y.polygon4land <- c(ucl4land, rev(lcl4land))

polygon(x.polygon4land, y.polygon4land, col = "#F2E399", border = NA)
lines(years4land, rowMeans(MTL4land), col = "#006349", lwd = 3)

fancyAxis(1, at = seq(1997, 2015, 2), las = 1, cex.axis = 2.5)
fancyAxis(2, at = seq(3.5, 4.7, 0.2), las = 1, cex.axis = 2.5)

rect(1998, 3.5, 2000, 4.7, col = "white", border = NA)


####################
# PLOT MTL_2.0 noP    #
####################

years2landnoP <- as.numeric(rownames(MTL2landnoP))
MTL2landnoP[is.na(MTL2landnoP)] <- 4

par(mar=c(6,7,2,2), mgp=c(4.5, 1, 0), las=0)

plot(NULL, xlim = c(min(years2landnoP), max(years2landnoP)), ylim=c(3.5,4.7),
     ylab="MTL_2.0",xlab="Year",cex.lab = 3,
     axes=FALSE, frame.plot = TRUE, bty = "l")

stdevs2landnoP <- apply(MTL2landnoP, 1, sd)
lcl2landnoP <- rowMeans(MTL2landnoP) - stdevs2landnoP
ucl2landnoP <- rowMeans(MTL2landnoP) + stdevs2landnoP
x.polygon2landnoP <- c(years2landnoP, sort(years2landnoP, decreasing = TRUE))
y.polygon2landnoP <- c(ucl2landnoP, rev(lcl2landnoP))

polygon(x.polygon2landnoP, y.polygon2landnoP, col = "#F2E399", border = NA)
lines(years2landnoP, rowMeans(MTL2landnoP), col = "#006349", lwd = 3)

fancyAxis(1, at = seq(1997, 2015, 2), las = 1, cex.axis = 2.5)
fancyAxis(2, at = seq(3.5, 4.7, 0.2), las = 1, cex.axis = 2.5)

rect(1998, 3.5, 2000, 4.7, col = "white", border = NA)


####################
# PLOT MTL_3.25 noP   #
####################

years3landnoP <- as.numeric(rownames(MTL3landnoP))
MTL3landnoP[is.na(MTL3landnoP)] <- 4

par(mar=c(6,7,2,2), mgp=c(4.5, 1, 0), las=0) 

plot(NULL, xlim = c(min(years3landnoP), max(years3landnoP)), ylim=c(3.5,4.7),
     ylab="MTL_3.25",xlab="Year",cex.lab = 3,
     axes=FALSE, frame.plot = TRUE, bty = "l")

stdevs3landnoP <- apply(MTL3landnoP, 1, sd)
lcl3landnoP <- rowMeans(MTL3landnoP) - stdevs3landnoP
ucl3landnoP <- rowMeans(MTL3landnoP) + stdevs3landnoP
x.polygon3landnoP <- c(years3landnoP, sort(years3landnoP, decreasing = TRUE))
y.polygon3landnoP <- c(ucl3landnoP, rev(lcl3landnoP))

polygon(x.polygon3landnoP, y.polygon3landnoP, col = "#F2E399", border = NA)
lines(years3landnoP, rowMeans(MTL3landnoP), col = "#006349", lwd = 3)

fancyAxis(1, at = seq(1997, 2015, 2), las = 1, cex.axis = 2.5)
fancyAxis(2, at = seq(3.5, 4.7, 0.2), las = 1, cex.axis = 2.5)

rect(1998, 3.5, 2000, 4.7, col = "white", border = NA)


####################
# PLOT MTL_4.0  noP  #
####################

years4landnoP <- as.numeric(rownames(MTL4landnoP))
MTL4landnoP[is.na(MTL4landnoP)] <- 4.2

par(mar=c(6,7,2,2), mgp=c(4.5, 1, 0), las=0) 

plot(NULL, xlim = c(min(years4landnoP), max(years4landnoP)), ylim=c(3.5,4.7),
     ylab="MTL_4.0",xlab="Year",cex.lab = 3,
     axes=FALSE, frame.plot = TRUE, bty = "l")

stdevs4landnoP <- apply(MTL4landnoP, 1, sd)
lcl4landnoP <- rowMeans(MTL4landnoP) - stdevs4landnoP
ucl4landnoP <- rowMeans(MTL4landnoP) + stdevs4landnoP
x.polygon4landnoP <- c(years4landnoP, sort(years4landnoP, decreasing = TRUE))
y.polygon4landnoP <- c(ucl4landnoP, rev(lcl4landnoP))

polygon(x.polygon4landnoP, y.polygon4landnoP, col = "#F2E399", border = NA)
lines(years4landnoP, rowMeans(MTL4landnoP), col = "#006349", lwd = 3)

fancyAxis(1, at = seq(1997, 2015, 2), las = 1, cex.axis = 2.5)
fancyAxis(2, at = seq(3.5, 4.7, 0.2), las = 1, cex.axis = 2.5)

rect(1998, 3.5, 2000, 4.7, col = "white", border = NA)


#########################################################################################
############################################################################################
###################  Statistical tests ####################################
############################################################################################
########################################################################################

##########################################################################
################ LM analysis (for parametric values) for landings ###########
##########################################################################

########MTL_2.0###########

MTL2land.t<-t(MTL2land)
MTL2land.t.m<-colMeans(MTL2land.t)
year<-c(1997:2014)
year<-as.numeric(year)
MTL2land.t.m<-as.data.frame(MTL2land.t.m)
MTL2land.t.m$year<-year
MTL_2.0b=as.data.frame(MTL2land.t.m)
View(MTL_2.0b)
MTL_2.0<-MTL_2.0b[-3,]
View(MTL_2.0)

plot(MTL_2.0$year, MTL_2.0$MTL2land.t.m, ylab="MTL_2.0", xlab="Years")
lm1<-lm(MTL_2.0$MTL2land.t.m ~ MTL_2.0$year, data=MTL_2.0)
shapiro.test(lm1$residuals) #normality if pvalue > alpha=0.05
library(lmtest)
hmctest(lm1) #if pvalue > alpha -> verified condition #Harrison-McCabe test statistic is the fraction of the residual sum of squares that relates to the fraction of the data before the breakpoint
dwtest(lm1) #no autocorrelation if p>alpha # checking independance (temporal autocorrelation)
summary(lm1) #significant slope if p<0.05

########MTL_3.25###########
MTL_3.25land.t<-t(MTL3land)
MTL_3.25land.t.m<-colMeans(MTL_3.25land.t)
year<-c(1997:2014)
year<-as.numeric(year)
MTL_3.25land.t.m<-as.data.frame(MTL_3.25land.t.m)
MTL_3.25land.t.m$year<-year
MTL_3.25b=as.data.frame(MTL_3.25land.t.m)
View(MTL_3.25b)
MTL_3.25<-MTL_3.25b[-3,]
View(MTL_3.25)

plot(MTL_3.25$year, MTL_3.25$MTL_3.25land.t.m, ylab="MTL_3.25", xlab="Years")
lm3<-lm(MTL_3.25$MTL_3.25land.t.m ~ MTL_3.25$year, data=MTL_3.25)
shapiro.test(lm3$residuals)
hmctest(lm3)
dwtest(lm3)
summary(lm3)


########MTL_4.0###########
MTL_4.0land.t<-t(MTL4land)
MTL_4.0land.t.m<-colMeans(MTL_4.0land.t)
year<-c(1997:2014)
year<-as.numeric(year)
MTL_4.0land.t.m<-as.data.frame(MTL_4.0land.t.m)
MTL_4.0land.t.m$year<-year
MTL_4.0b=as.data.frame(MTL_4.0land.t.m)
View(MTL_4.0b)
MTL_4.0<-MTL_4.0b[-3,]
View(MTL_4.0)

plot(MTL_4.0$year, MTL_4.0$MTL_4.0land.t.m, ylab="MTL_4.0", xlab="Years")
lm9<-lm(MTL_4.0$MTL_4.0land.t.m ~ MTL_4.0$year, data=MTL_4.0)
shapiro.test(lm9$residuals) 
hmctest(lm9)
dwtest(lm9)
summary(lm9)

########MTL_2.0 noP###########

MTL2landnoP.t<-t(MTL2landnoP)
MTL2landnoP.t.m<-colMeans(MTL2landnoP.t)
year<-c(1997:2014)
year<-as.numeric(year)
MTL2landnoP.t.m<-as.data.frame(MTL2landnoP.t.m)
MTL2landnoP.t.m$year<-year
MTL_2.0c=as.data.frame(MTL2landnoP.t.m)
View(MTL_2.0c)
MTL_2.0noP<-MTL_2.0c[-3,]
View(MTL_2.0noP)

plot(MTL_2.0noP$year, MTL_2.0noP$MTL2landnoP.t.m, ylab="MTL_2.0", xlab="Years")
lm2<-lm(MTL_2.0noP$MTL2landnoP.t.m ~ MTL_2.0noP$year, data=MTL_2.0noP)
shapiro.test(lm2$residuals)
library(lmtest)
hmctest(lm2)
dwtest(lm2) 
summary(lm2) 

########MTL_3.25 noP###########
MTL_3.25landnoP.t<-t(MTL3landnoP)
MTL_3.25landnoP.t.m<-colMeans(MTL_3.25landnoP.t)
year<-c(1997:2014)
year<-as.numeric(year)
MTL_3.25landnoP.t.m<-as.data.frame(MTL_3.25landnoP.t.m)
MTL_3.25landnoP.t.m$year<-year
MTL_3.25c=as.data.frame(MTL_3.25landnoP.t.m)
View(MTL_3.25c)
MTL_3.25noP<-MTL_3.25c[-3,]
View(MTL_3.25noP)

plot(MTL_3.25noP$year, MTL_3.25noP$MTL_3.25landnoP.t.m, ylab="MTL_3.25", xlab="Years")
lm4<-lm(MTL_3.25noP$MTL_3.25landnoP.t.m ~ MTL_3.25noP$year, data=MTL_3.25noP)
shapiro.test(lm4$residuals)
hmctest(lm4)
dwtest(lm4)
summary(lm4)


########MTL_4.0 noP###########
MTL_4.0landnoP.t<-t(MTL4landnoP)
MTL_4.0landnoP.t.m<-colMeans(MTL_4.0landnoP.t)
year<-c(1997:2014)
year<-as.numeric(year)
MTL_4.0landnoP.t.m<-as.data.frame(MTL_4.0landnoP.t.m)
MTL_4.0landnoP.t.m$year<-year
MTL_4.0c=as.data.frame(MTL_4.0landnoP.t.m)
View(MTL_4.0c)
MTL_4.0noP<-MTL_4.0c[-3,]
View(MTL_4.0noP)

plot(MTL_4.0noP$year, MTL_4.0noP$MTL_4.0landnoP.t.m, ylab="MTL_4.0", xlab="Years")
lm6<-lm(MTL_4.0noP$MTL_4.0landnoP.t.m ~ MTL_4.0noP$year, data=MTL_4.0noP)
shapiro.test(lm6$residuals) 
hmctest(lm6)
dwtest(lm6)
summary(lm6)

#############################################################################
################ MannKendall analysis for landings ###########
##########################################################################

##MTL_2.0#####
acf(MTL_2.0)
pacf(MTL_2.0)
library(Kendall)
res2 <-  MannKendall(MTL_2.0$MTL2land.t.m)
summary(res2)
#p-value associated with the Mann-Kendall test is statistically significant when p<0.05. 
#Tau  -> Kendall's tau statistic
#p-value   -> two-sided p-value
#Score	-> Kendall Score (S)
#Denominator	-> denominator (D), such that tau=S/D
#Var(Score) 	-> variance of Kendall Score

##MTL_3.25#####
acf(MTL_3.25)
pacf(MTL_3.25)
res3 <-  MannKendall(MTL_3.25$MTL_3.25land.t.m)
summary(res3)

##MTL_4.0#####
acf(MTL_4.0)
pacf(MTL_4.0)
res4 <-  MannKendall(MTL_4.0$MTL_4.0land.t.m)
summary(res4)


##MTL_2.0 noP#####
acf(MTL_2.0noP)
pacf(MTL_2.0noP)
res5 <-  MannKendall(MTL_2.0noP$MTL2landnoP.t.m)
summary(res5)

##MTL_3.25 noP#####
acf(MTL_3.25noP)
pacf(MTL_3.25noP)
res6 <-  MannKendall(MTL_3.25noP$MTL_3.25landnoP.t.m)
summary(res6)

##MTL_4.0 noP#####
acf(MTL_4.0noP)
pacf(MTL_4.0noP)
res7 <-  MannKendall(MTL_4.0noP$MTL_4.0landnoP.t.m)
summary(res7)
