# This script (1) allows the update of the classification of trends based on the new thresholds
# This script (2) allows the creation of an alluvial plot based on switch in trend classification
# It is based on changes on 2 functions : classifyLinearTrends and plotLinearClassification
# It updates : 
# - the figures obtained showing all species classification trends (points with error bar or barplot)
# - the tables obtained especially the column "trend"
# NB : the script works if the working directory is the root of the indicator routine

rm(list = ls()) 

# Load libraries
library(tidyverse)
library(ggalluvial)

# Result path
result_path = "C:/Users/mvimont01/Downloads/maj_appendix_plants/"

# FUNCTIONS ----
# Load all functions required from the routine
source("R/identifyConvIssue.R")
source("R/formatTrendEstimates.R")
source("R/plotLinearClassification.R")
source("R/fillAbsence.R")

# Updated version of the function 'classifyLinearTrends' 
classifyLinearTrends_bis = function(data, distribution, threshold_inf = 0.964, threshold_sup = 1.036){
  data$trend = NA
  
  # Erase NA values
  dataNA = data[is.na(data$estimate),]
  data = data[!is.na(data$estimate),]
  
  # Turn to exponential estimates and confidence intervals
  if(distribution != "gaussian"){
    data[,c("estimate", "infIC", "supIC")] = exp(data[,c("estimate", "infIC", "supIC")])
    limit = 1
  }else{
    limit = 0
  }
  
  # Decrease
  data$trend[data$pval < 0.05 & data$estimate < limit] = "Déclin modéré"
  data$trend[data$pval < 0.05 & data$estimate < threshold_inf] = "Déclin modéré à fort"
  data$trend[data$pval < 0.05 & data$supIC < threshold_inf] = "Déclin fort"
  
  # Increase
  data$trend[data$pval < 0.05 & data$estimate > limit] = "Augmentation modérée"
  data$trend[data$pval < 0.05 & data$estimate > threshold_sup] = "Augmentation modérée à forte"
  data$trend[data$pval < 0.05 & data$infIC > threshold_sup] = "Augmentation forte"
  
  # Stability
  data$trend[data$pval > 0.05 & (data$infIC < threshold_inf | data$supIC > threshold_sup)] = "Incertain"
  data$trend[data$pval > 0.05 & data$infIC >= threshold_inf & data$supIC <= threshold_sup] = "Stable"
  
  # Turn back coefficients
  if(distribution != "gaussian"){
    data[,c("estimate", "infIC", "supIC")] = log(data[,c("estimate", "infIC", "supIC")])
  }
  
  # Bind all dataframes
  data = rbind(data, dataNA)
  
  return(data)
}


# Updated version of the function 'plotLinearClassification' 
plotLinearClassification_bis <- function(data, distribution, type = "barplot", threshold_inf = 0.964, threshold_sup = 1.036,
                                         path = ""){
  
  # Erase NA values
  data = data[!is.na(data$trend),]
  
  # Transform coefficients according to distribution
  if(distribution != "gaussian"){
    data[,c("estimate", "infIC", "supIC")] = exp(data[,c("estimate", "infIC", "supIC")])
    zeroVal = 1
  }else{
    zeroVal = 0
  }
  
  # Order and colors for trends
  ## Names
  refOrder = c("Strong decline", "Moderate to strong decline", "Moderate decline", 
               "Stable", "Moderate increase", "Moderate to strong increase", 
               "Strong increase", "Uncertain")
  
  ## Colors
  refCols = c(rev(c("#274034", "#4d8067", "#81b29a", "#F2CC8F","#f57e56","#d8400d", "#d61a0d")),  "gray")
  
  ## Point types
  refTypes = c(16, 1, 3, 3, 3, 1, 16, 4)
  
  # Filter for trends that occur
  refCols = refCols[which(!is.na(match(refOrder,unique(data$trend))))]
  refTypes = refTypes[which(!is.na(match(refOrder,unique(data$trend))))]
  refOrder = refOrder[which(!is.na(match(refOrder,unique(data$trend))))]
  
  # Reorder data trends with refOrder
  data$trend = factor(data$trend, levels = refOrder)
  
  # Reorder data species with decreasing growthrate
  data = dplyr::arrange(data, estimate)
  data$species = factor(data$species, data$species)
  
  # Make barplot
  if(type == "barplot"){
    plot = ggplot2::ggplot(data, ggplot2::aes(x = trend, fill = trend)) +
      
      # Make barplot
      ggplot2::geom_bar(stat = "count") +
      
      # Adapt colors
      ggplot2::stat_count(geom = "text", colour = "gray20", size = 6.5,
                          ggplot2::aes(y = ggplot2::after_stat(count) + 8, label = ggplot2::after_stat(count))) +
      
      # Change colors
      ggplot2::scale_fill_manual("Classification", values = refCols) +
      
      # Change theme
      ggplot2::theme_bw() + 
      
      # Set a maximum y value
      ggplot2::ylim(c(0,max(table(data$trend))+10)) + 
      
      # Change axis names
      ggplot2::xlab("") + ggplot2::ylab("Number of species") +
      
      # Format theme
      ggplot2::theme(panel.grid.major.x = ggplot2::element_blank(),
                     panel.grid.minor.y = ggplot2::element_blank(),
                     axis.text.x = ggplot2::element_blank(), 
                     axis.ticks.x = ggplot2::element_blank(),
                     axis.text.y = ggplot2::element_text(size = 18),
                     axis.title.y = ggplot2::element_text(size = 20),
                     legend.text = ggplot2::element_text(size = 12),
                     legend.title = ggplot2::element_text(size = 14),
                     legend.position = "none")
    
  }else if(type == "errorbar"){
    plot = ggplot2::ggplot(data, ggplot2::aes(y = species, x = estimate, col = trend, shape = trend)) + 
      # Add estimates and errorbars
      ggplot2::geom_pointrange(ggplot2::aes(xmin = infIC, xmax = supIC), linewidth = 0.8, fatten = 3) +
      
      # Add vertical line for stable/uncertain effect
      ggplot2::geom_vline(xintercept = zeroVal, linetype = "dashed", col = "gray") +
      
      # Add vertical line for strong negative trend
      ggplot2::geom_vline(xintercept = threshold_inf, linetype = "dashed", col = "#d61a0d") +
      
      # Add vertical line for strong positive effect
      ggplot2::geom_vline(xintercept = threshold_sup, linetype = "dashed", col = "#274034") +
      
      # Change color
      ggplot2::scale_color_manual("Classification", values = refCols) +
      
      # Change point types
      ggplot2::scale_shape_manual("Classification", values = refTypes) +
      
      # Change names of axis
      ggplot2::ylab("Species") + ggplot2::xlab("Growth rates") + 
      
      # Apply classic theme
      ggplot2::theme_bw() +
      
      # Format theme
      ggplot2::theme(axis.ticks.y = ggplot2::element_blank(), 
                     axis.text.y = ggplot2::element_blank(),
                     axis.title = ggplot2::element_text(size = 20),
                     legend.title = ggplot2::element_text(size = 14),
                     legend.text = ggplot2::element_text(size = 12),
                     axis.text.x = ggplot2::element_text(size = 18),
                     legend.position = "none")
    
  }
  
  # Save result
  ggplot2::ggsave(plot = plot, filename = paste0(type, ".png"), device = "png", path = path, 
                  width = 20, height = 15, units = "cm")
  
  return(plot)
}

# 1. UPDATING THE FIGURES  ----
# Import the presence/absence data.Frame 
dataMT = data.table::fread(file = "data/VF_2023/countingData.csv", encoding = "UTF-8", drop = "V1")
dataMT = data.frame(dataMT[,c("ID","species","point", "year", "site", "julian", "r_resp", "abond_cumul", "fail")])

# Extract the list of species analyzed in the paper
speciesList = c("Hedera helix", "Lolium perenne", "Rubus fruticosus",
                 "Dactylis glomerata", "Galium aparine", "Plantago lanceolata",
                 "Convolvulus arvensis", "Quercus robur", "Fraxinus excelsior",
                 "Trifolium repens", "Polygonum aviculare", "Fallopia convolvulus",
                 "Carpinus betulus", "Crataegus monogyna", "Holcus lanatus",
                 "Urtica dioica", "Arrhenatherum elatius", "Cirsium arvense",
                 "Geranium dissectum", "Daucus carota", "Geum urbanum",
                 "Lysimachia arvensis", "Lonicera periclymenum", "Chenopodium album",
                 "Poa annua",  "Prunus spinosa", "Geranium robertianum",
                 "Plantago major", "Picris hieracioides", "Anisantha sterilis",
                 "Achillea millefolium", "Trifolium pratense", "Acer pseudoplatanus",
                 "Poa trivialis", "Veronica persica", "Heracleum sphondylium",
                 "Potentilla reptans", "Medicago lupulina", "Fagus sylvatica",
                 "Pteridium aquilinum", "Cornus sanguinea", "Corylus avellana",
                 "Ranunculus repens", "Vicia sativa", "Rubia peregrina",
                 "Ligustrum vulgare", "Sonchus asper", "Brachypodium sylvaticum",
                 "Glechoma hederacea", "Alopecurus myosuroides","Acer campestre",
                 "Clematis vitalba", "Poa pratensis", "Hypochaeris radicata",
                 "Lotus corniculatus", "Viola arvensis", "Avena fatua",
                 "Geranium molle", "Castanea sativa", "Rumex acetosa",
                 "Hypericum perforatum", "Bellis perennis", "Papaver rhoeas",
                 "Taraxacum officinale", "Senecio vulgaris", "Bromus hordeaceus",
                 "Elytrigia repens", "Lapsana communis", "Convolvulus sepium",
                 "Ranunculus acris", "Veronica arvensis", "Agrostis stolonifera",
                 "Tripleurospermum inodorum","Solanum nigrum", "Sonchus oleraceus",
                 "Carex sylvatica", "Agrostis capillaris", "Fragaria vesca",
                 "Equisetum arvense", "Rubus caesius", "Cerastium fontanum",
                 "Prunella vulgaris", "Teucrium scorodonia", "Stellaria media",
                 "Galium mollugo", "Poa nemoralis", "Poterium sanguisorba",
                 "Quercus ilex", "Ulmus minor",  "Quercus petraea",
                 "Jacobaea vulgaris", "Ilex aquifolium", "Prunus avium",
                 "Anthoxanthum odoratum", "Lolium multiflorum", "Persicaria maculosa",
                 "Artemisia vulgaris", "Quercus pubescens", "Euonymus europaeus",
                 "Erigeron canadensis", "Festuca rubra", "Cirsium vulgare",
                 "Veronica chamaedrys", "Agrimonia eupatoria", "Crepis capillaris",
                 "Matricaria chamomilla", "Capsella bursa-pastoris","Polygonatum multiflorum",
                 "Rumex obtusifolius", "Rosa canina",  "Schedonorus arundinaceus",
                 "Silene latifolia", "Viola reichenbachiana", "Ervilia hirsuta",
                 "Geranium rotundifolium","Helminthotheca echioides", "Molinia caerulea",
                 "Ruscus aculeatus", "Rumex crispus", "Lactuca serriola",
                 "Myosotis arvensis", "Geranium pusillum", "Holcus mollis",
                 "Asparagus acutifolius", "Origanum vulgare", "Stellaria holostea",
                 "Carex flacca", "Alliaria petiolata", "Calluna vulgaris",
                 "Mercurialis annua", "Cytisus scoparius", "Leucanthemum vulgare",
                 "Phillyrea angustifolia","Pinus sylvestris", "Mercurialis perennis",
                 "Triticum aestivum", "Euphorbia amygdaloides","Melica uniflora",
                 "Sherardia arvensis", "Dioscorea communis", "Ranunculus bulbosus",
                 "Teucrium chamaedrys", "Hordeum murinum", "Anemone nemorosa",
                 "Anthriscus sylvestris", "Smilax aspera", "Epilobium tetragonum",
                 "Thymus vulgaris", "Buxus sempervirens", "Medicago sativa",
                 "Circaea lutetiana", "Frangula alnus", "Phleum pratense",
                 "Rubus ulmifolius", "Robinia pseudoacacia", "Ajuga reptans",
                 "Vicia sepium", "Acer platanoides", "Euphorbia cyparissias",
                 "Brachypodium retusum", "Aphyllanthes monspeliensis", "Betula pendula",
                 "Bromopsis erecta", "Sinapis arvensis", "Viburnum lantana",
                 "Brachypodium pinnatum", "Sambucus nigra", "Pilosella officinarum",
                 "Potentilla erecta", "Pastinaca sativa", "Ulex europaeus",
                 "Convallaria majalis", "Sorbus torminalis", "Viola odorata",
                 "Galeopsis tetrahit", "Hyacinthoides non-scripta", "Moehringia trinervia",
                 "Rosa arvensis", "Zea mays", "Quercus coccifera",
                 "Sorbus aucuparia")


# Add 0's (absences) to the data.frame
dataMT = fillAbsence(dataMT, interestVar = c("abond_cumul", "fail"), speciesList, method = "once")

# Extract data for the 10 latest years
dataST = dataMT[dataMT$year >= 2013,]

# Format trend estimations for both long and short-term trends
dataMediumTerm = formatTrendEstimates(dataMT, speciesList, "VF_2023", "longTermTrend")
dataShortTerm = formatTrendEstimates(dataST, speciesList, "VF_2023", "shortTermTrend")

# Classify those trends into 8 classes
medClass = classifyLinearTrends_bis(data = dataMediumTerm, distribution = "betabinomial", threshold_inf = 0.961, threshold_sup = 1.030)
shortClass = classifyLinearTrends_bis(data = dataShortTerm, distribution = "betabinomial", threshold_inf = 0.961, threshold_sup = 1.030)

# Make the figures fot the long-term trends
medClass$trend = factor(medClass$trend, levels  = c("Augmentation forte","Augmentation modérée à forte", "Augmentation modérée", "Stable",
                                                    "Déclin modéré", "Déclin modéré à fort", "Déclin fort", "Incertain"),
                        labels =  c("Strong increase", "Moderate to strong increase", "Moderate increase",
                                    "Stable", "Moderate decline", "Moderate to strong decline", "Strong decline", 
                                    "Uncertain"))
plotLinearClassification_bis(medClass, "betabinomial", type = "barplot", threshold_inf = 0.961, threshold_sup = 1.030, path = "C:/Users/mvimont01/Downloads/maj_appendix_plants/")
plotLinearClassification_bis(medClass, "betabinomial", type = "errorbar", threshold_inf = 0.961, threshold_sup = 1.030, path = "C:/Users/mvimont01/Downloads/maj_appendix_plants/")

# Make the figures fot the short-term trends
shortClass$trend = factor(shortClass$trend, levels  = c("Augmentation forte","Augmentation modérée à forte", "Stable",
                                                        "Déclin modéré à fort", "Déclin fort", "Incertain"),
                        labels =  c("Strong increase", "Moderate to strong increase", 
                                    "Stable", "Moderate to strong decline", "Strong decline", 
                                    "Uncertain"))
plotLinearClassification_bis(shortClass, "betabinomial", type = "barplot",  threshold_inf = 0.961, threshold_sup = 1.030, path = "C:/Users/mvimont01/Downloads/maj_appendix_plants/")
plotLinearClassification_bis(shortClass, "betabinomial", type = "errorbar",  threshold_inf = 0.961, threshold_sup = 1.030, path = "C:/Users/mvimont01/Downloads/maj_appendix_plants/")



# 2. UPDATING THE TABLES  ----

# Extract trends based on the new classification thresholds 
medClass = medClass[,c("species", "trend")]
shortClass = shortClass[,c("species", "trend")]

# Import old trends
dataTrend_MT = read.csv("outputs/VF_2023/tables/trends/longTermTrends.csv", row.names = 1)
dataTrend_ST = read.csv("outputs/VF_2023/tables/trends/shortTermTrends.csv", row.names = 1)

# Erase columns that are no more relevant
dataTrend_MT = dataTrend_MT[,-c(16:18)] 
dataTrend_ST = dataTrend_ST[,-c(16:18)] 

# Make and save the new table (with updated trend classification) 
dataTrend_MT = merge(dataTrend_MT, medClass, by = "species")
write.csv(dataTrend_MT, paste0(result_path,"longTermTrends_V2.csv"))
dataTrend_ST = merge(dataTrend_ST, shortClass, by = "species")
write.csv(dataTrend_ST, paste0(result_path, "shortTermTrends_V2.csv"))

# 3. ALLUVIAL PLOT ----
# Adapt colnames for alluvial
colnames(medClass)[2] = "trend_MT"
colnames(shortClass)[2] = "trend_ST"

# Merge both long and short-term trends
dataClass = merge(medClass, shortClass, by = "species")

# Create the alluvial dataframe by counting changes in trends
dataAlluvial = group_by(dataClass, trend_MT, trend_ST) %>%
  summarise(count = n())

# Erase trends that cannot be calculated
dataAlluvial = dataAlluvial[!is.na(dataAlluvial$trend_MT) & !is.na(dataAlluvial$trend_ST),]

# Erase trends that are uncertain in both temporal frames
dataAlluvial = dataAlluvial[!(dataAlluvial$trend_MT == "Uncertain" & dataAlluvial$trend_ST == "Uncertain"),]

# Reorder the classification scheme
lab = c("Strong increase", "Moderate to strong increase", "Moderate increase", 
        "Stable", "Moderate decline", "Moderate to strong decline", "Strong decline", 
        "Uncertain")
dataAlluvial$trend_ST = factor(dataAlluvial$trend_ST, levels = levels(dataAlluvial$trend_MT)[-c(3,5)], labels = lab[-c(3,5)])
dataAlluvial$trend_MT = factor(dataAlluvial$trend_MT, levels = levels(dataAlluvial$trend_MT), labels = lab)

# Chose colors
colTrends = c("#274034", "#4d8067", "#81b29a", "#F2CC8F","#f57e56","#d8400d", "#d61a0d",  "gray")

# Make plot
ggplot(dataAlluvial,
       aes(y = count, axis1 = trend_MT, axis2 = trend_ST)) +
  geom_alluvium(data = dataAlluvial, aes(fill = trend_MT), alpha = 0.9) + 
  geom_stratum(fill = rev(c(colTrends[-c(3,5)],colTrends)),
               width = 1/6, col = "gray10") +
  scale_fill_manual("Classification", values = colTrends) +
  theme_bw() +
  geom_text(stat = "stratum", aes(label = after_stat(stratum)),
            reverse = TRUE, size = 2.5, col = "black") +
  ylab("Number of species") +
  scale_x_discrete(limits = c("2009 - 2023\nperiod", "2013 - 2023\nperiod"), expand = c(.05, .05)) +
  theme(legend.position = "bottom",
        axis.text.x = element_text(angle = 45, hjust = 1, size = 14),
        axis.title.y = element_text(size = 16),
        axis.text.y = element_text(size = 14))

# Save plot
ggsave("alluvial_plot_zoom_VF_V2.png", last_plot(), path = result_path, device = "png", 
       width = 22, height = 18, units = "cm")
