1_0_Create SS_Summary_and_Evaluation.Rmd

---
title: "Check sus and prepare summarized vegetation data by site units"
author: "Will MacKenzie"
date: "13/09/2021"
output: html_document
---

```{r setup, include=FALSE}
require(data.table)
require(tidyverse)
require(dplyr)
require(data.tree)
require(DBI) #loads odbc as well
require(labdsv)
require(factoextra)
require(CooccurrenceAffinity)
require(dendextend)
```

#### Import raw data and convert into analysis set
```{r set folders, include=FALSE}
#source("./_functions/_convert_vpro_veg.R")
#  becmaster <- dbConnect(odbc::odbc(), .connection_string = "Driver={Microsoft Access Driver (*.mdb, *.accdb)}; DBQ=D:/GitHub/BECMaster_Cleaning/updated_vpro/BECMaster_fixing.accdb;")
# 
# plot.veg <- dbReadTable(becmaster, "BECMaster_fixing_Veg")
# 
# dbDisconnect(becmaster)
# 
# sppmaster <- dbConnect(odbc::odbc(), .connection_string = "Driver={Microsoft Access Driver (*.mdb, *.accdb)}; DBQ=F:/OneDrive - Personal/OneDrive/BCSpeciesList/SpeciesTaxonomyMaster.accdb;")
# taxon.all  <- dbReadTable(sppmaster, "USysAllSpecs")
# dbDisconnect(sppmaster)
# taxon.lifeform <- taxon.all %>% filter(Codetype == "U" |Codetype == "X") %>% dplyr::select(Code, Lifeform) %>% distinct
# veg.dat <- convert_vpro_veg(plot.veg, taxon.lifeform)
# save(veg.dat, file = "./clean_data/Analysis_BECMaster_Veg.RData")
```

#### Import all analysis data and other vpro tables for building summary
```{r set folders, include=FALSE}
veg.dat <- readRDS("./clean_data/Analysis_BECMaster_Veg.rds") ###named veg.dat

 becmaster <- dbConnect(odbc::odbc(), .connection_string = "Driver={Microsoft Access Driver (*.mdb, *.accdb)}; DBQ=D:/GitHub/BECMaster_Cleaning/updated_vpro/BECMaster_fixing.accdb;")
plot.env <- dbReadTable(becmaster, "BECMaster_fixing_Env")
dbDisconnect(becmaster)

# ###created in 0_0_Plot_su_Hierarchy checks script
master_su <- dbConnect(odbc::odbc(), .connection_string = "Driver={Microsoft Access Driver (*.mdb, *.accdb)};
DBQ=D:/BC_Correlation2_Vpro_2023/CoastGuide_Hierarchy.accdb;")
su <- dbReadTable(master_su, "All_Coast_Forest_2024v4_SU")
hier <- dbReadTable(master_su, "CoastForest_v2024_2_Hierarchy")
dbDisconnect(master_su)

# master_su <- dbConnect(odbc::odbc(), .connection_string = "Driver={Microsoft Access Driver (*.mdb, *.accdb)};
# DBQ=D:/BC_Correlation2_Vpro_2023/working/All_Skunkcabbage.accdb;")
# su <- dbReadTable(master_su, "AllBGC_Skunkcabbage_SU")
# hier <- dbReadTable(master_su, "AllBGC_Skunkcabbage_Hierarchy")
# dbDisconnect(master_su)


sppmaster <- dbConnect(odbc::odbc(), .connection_string = "Driver={Microsoft Access Driver (*.mdb, *.accdb)}; DBQ=F:/OneDrive - Personal/OneDrive/BCSpeciesList/SpeciesTaxonomyMaster.accdb;")
taxon.all  <- dbReadTable(sppmaster, "USysAllSpecs")
dbDisconnect(sppmaster)
taxon.lifeform <- taxon.all %>% filter(Codetype == "U" |Codetype == "X") %>% dplyr::select(Code, Lifeform) %>% distinct

veglump <- dbConnect(odbc::odbc(), .connection_string = "Driver={Microsoft Access Driver (*.mdb, *.accdb)}; DBQ=D:/BC_Correlation2_Vpro_2023/CoastGuide_Spp_lump.accdb;")
lump <- dbReadTable(veglump, "CoastGuide2023_lump")
dbDisconnect(veglump) 
```


## Import su table and BECdb for current BGC units and site series. Look for missing units. Create complete SU
Working on Coast guide first

```{r import su tables and look for errors}
###remove phases and seral units
su_siteseries <- su %>% 
  filter(!str_detect(SiteUnit, '[$]'))
phases = c("a", "b", "c")
su_siteseries$SiteUnit <- str_replace(su_siteseries$SiteUnit, "[abc]", "")
fwrite(su_siteseries, "./clean_tabs/CoastGuide_All_BGC_Forest_SS.csv")

su <- su_siteseries
###Check that there is plot data for all
missingplots <- anti_join(su, plot.env)
plots <- right_join(su, plot.env)
### Check for duplicate plot usage
usemorethanone <- su %>% dplyr::group_by(PlotNumber) %>% dplyr::mutate(dups = n()) %>% filter(dups >1) %>% ungroup() %>% arrange(PlotNumber)
fwrite(usemorethanone, "./review_outputs/Coast_Plots_Used_Morethan_Once.csv")
###Summaries by BGC and by SS
#plots_bgc <- su %>% group_by(bgc) %>% dplyr::summarise(plots = n())
plots_SS <- su %>% group_by(SiteUnit) %>% dplyr::summarise(plots = n())
SS_count <- length(unique(su$SiteUnit))
siteunits_toofew <- plots_SS %>% filter(plots <5) #%>% dplyr::select(- PlotNumber) %>% distinct

fwrite(siteunits_toofew, "./review_outputs/Coast_SiteUnits_w_lessthan5plots.csv")
```
#Roll up into hierarchy units and fopr running
```{r}
source("./_functions/_TabletoTree.R")
source("./_functions/_TreetoTable.R")
hier2 <- treeToTable(hier)
hier2 <- hier2[[1]]
hier3 <- hier2 %>% filter(Order == "Order Hw") %>% left_join(su)
su2 <- hier3 %>% select(PlotNumber, Assoc) %>% rename(SiteUnit = 2)
```


##Internal Metrics of Homogeneity of site series - Vegetation
Number of plots
Similarity between plots
Outliers

```{r internal variability}
require(tabula)
require(labdsv)
require(proxyC)
veg.dat <- veg.dat[PlotNumber %in% su$PlotNumber] %>% left_join(su)## recuce to plots in SU only
veg.unit <- veg.dat %>% filter(SiteUnit == "CWH vh 2 /101") %>% select(-Lifeform, -SiteUnit) %>% filter(Cover > .99) %>% matrify() %>% as.matrix 
veg.unit[veg.unit==0]<-NA
xx <- simil(veg.unit, method = "ejaccard", min_simil = .1) %>% as.matrix
xx

```

##roll up into site series summary data and reduce to analysis set
```{r summarize site series, echo=FALSE}
source("./_functions/_lump_species.R")
### Split into 3 indicator groups, tree, moss/lichen, and others
moss = c("9", "10", "11"); tree = c("1", "2"); understorey = c("3", "4", "5", "6", "7", "8","12")
vegdata.tree <- veg.dat %>% filter(Lifeform %in% tree) %>% lump_species(lump) %>% group_by(PlotNumber, Species) %>% mutate(totalcov = sum(Cover))
vegdata.understorey <- veg.dat %>% filter(Lifeform %in% understorey)%>% lump_species(lump)%>% group_by(PlotNumber, Species) %>% mutate(totalcov = sum(Cover))
vegdata.moss <- veg.dat %>% filter(Lifeform %in% moss) %>% lump_species(lump)%>% group_by(PlotNumber, Species) %>% mutate(totalcov = sum(Cover))

#vegdata2 <- lump_species(lump)

#vegdata <- vegdata %>% filter(lifeform <3)
source("./_functions/_create_su_vegdata.R")
source("./_functions/_create_analysis_vegsum.R")

veg_anal.tree <- create_su_vegdata(vegdata.tree, su2) %>% create_analysis_vegsum (minimportance = 0.5, minconstancy = 60, minplots = 5)#, covadj = 1)
veg_anal.under <- create_su_vegdata(vegdata.understorey, su2) %>% create_analysis_vegsum (minimportance = 0.5, minconstancy = 60, minplots = 5)#, covadj = .5)
veg_anal.moss <- create_su_vegdata(vegdata.moss, su2) %>% create_analysis_vegsum (minimportance = 0.5, minconstancy = 60, minplots = 5)#, covadj = .1)
veg_anal <- rbind(veg_anal.tree, veg_anal.under, veg_anal.moss)
# save(veg_anal, file = "./clean_data/SS_sum_analysis_all.RData")
# 
# save(veg_anal, file = "./clean_data/SS_sum_analysis_reduced.RData")
# SS_reduced_list <- unique(veg_anal$SiteUnit) %>% as.data.frame %>% dplyr::rename(SiteUnit = 1)
# su_reduced <- su %>% filter(SiteUnit %in% SS_reduced_list$SS)
```

Cluster of site units to create associations.
Uses the reduced species set and sets an  distance cut-off that seems to equate to association
Need to test against the pairwise routine.

```{r create distance matrix}
#require(harrietr)
#load("./clean_data/SS_sum_analysis_reduced.RData")
# test.su = c("CWH vm 1 /116.1", "CWH vm 3 /116.1")
# test.dat <- fread("./inputs/test.csv") %>% matrify
veg.dat <-
  veg_anal %>% dplyr::select(SiteUnit, Species, spp_importance) %>% filter(SiteUnit %in% test.su) %>% mutate(spp_importance = log10(spp_importance +1)) %>% matrify
 test_jac.vegan <- vegan::vegdist(test.dat, method = "bray", binary = F) %>% as.matrix(labels = TRUE)
test_jac.phil <- philentropy::distance(test.dat, method = "czekanowski", use.row.names = TRUE) %>% as.matrix(labels = TRUE)
test_jac <- proxy::dist(veg.dat, method="ejaccard", diag=FALSE, upper = FALSE, pairwise = FALSE) %>% as.matrix(labels = TRUE)## ejaccard and philentropy jaccard are real-data versions extended jaccard opr taimoto
xx <- as.matrix(test_jac)
# units <- as.data.frame(rownames(xx))

# , use.row.names = TRUE, as.dist.obj = TRUE)
#test_jac.tab <- data.frame(t(apply(test_jac, 1, sort))) %>% unique
test_jac.tab <- test_jac %>% as.matrix %>% as.data.frame
test_jac.tab[lower.tri(test_jac.tab)] <- 0
test_jac.list <- test_jac.tab %>%  rownames_to_column("SiteUnit1") %>% pivot_longer(!SiteUnit1) %>% rename(dist = value, SiteUnit2 = name) %>% as.data.frame %>% mutate_if(is.factor, as.character)# %>% filter(dist<.25, dist>0)

minimum.iso1 <- test_jac.list %>% group_by(SiteUnit1) %>%  slice_min(order_by = dist) %>% filter(dist<.25) %>% dplyr::select(SiteUnit1) %>% dplyr::rename(SS = SiteUnit1)
minimum.iso2 <- test_jac.list %>% group_by(SiteUnit2) %>%  slice_min(order_by = dist) %>% filter(dist<.25) %>% dplyr::select(SiteUnit2) %>% dplyr::rename(SS = SiteUnit2)
SS_w_neighbours <- rbind(minimum.iso1,  minimum.iso2)  %>% ungroup %>% unique %>% rename(SiteUnit = SS)
Root_su <- left_join(SS_w_neighbours, su) %>% select(PlotNumber, SiteUnit)
fwrite(Root_su , "./outputs/SS_to_merge.csv")
#SS_w_neighbours2 <- left_join(SS_w_neighbours, test_jac.list, by = c("SS" = "col")) %>% unique 
SS_no_neighbours <- anti_join(SS_reduced_list, SS_w_neighbours)

fwrite(test_jac.tab, "./outputs/dissim_matrix2.csv", row.names = TRUE)
fwrite(test_jac.list, "./outputs/dissim_pairwise2.csv")
fwrite(SS_no_neighbours, "./outputs/SS_with_no_similar.csv")
fwrite(SS_w_neighbours , "./outputs/SS_with_no_similar.csv")
```

Use agnes algorithms to build bottom-up hierarchy. Cutting the dendrogram at a low dissimilarity threshold (0.25) seems to make good working groups that pass the second pair-wise test. However, the structure of the higher level groups seems odd in some cases where similar site series end up in very different upper units.
Might need to develop an alternative grouping mechanism
```{r cluster and cut - standard algorithms}
#ss_clst <- hclust(test_jac, method = "complete")
library(cluster)
require(ape)
ss_clst = agnes(test_jac, diss = TRUE, stand = FALSE,
      method = "complete")
# ht_dendro <- max(ss_clst$height)*.25
# dendro_test <- ss_clst %>% as.dendrogram %>%
#      set("branches_k_color", k=8) %>% set("branches_lwd", c(1)) %>%
#    set("branches_lty", c(1)) %>%
#    set("labels_colors") %>% set("labels_cex", c(.5)) %>% 
#    set("nodes_pch", 19) %>% set("nodes_col", c("black"))
# plot(dendro_test)

dendro_hc <- as.hclust(ss_clst)
fviz_dend(dendro_hc, cex = 0.5, lwd = 0.5, h = .25, 
          rect = TRUE, 
          k_colors = "jco", 
          rect_border = "jco", 
          rect_fill = TRUE,
          ggtheme = theme_gray(),labels=F)

# ggplot(dendro_test, horiz = TRUE); abline(v = .25, col = 2, lty = 2)
# 
# plot(horiz = TRUE)
# 
# 
# 
# cutgrps <- cutree(dendro_hc, h=.25)
# plot(as.phylo(dendro_hc), cex = .5, label.offset = .01)
# 
# par(mfrow=c(3,1))
# 
# plot(dendro_test, main="Main")
# plot(cut(dendro_test, h = ht_dendro)$upper,
#      main="Upper tree of cut at h=25")
# plot(cut(dendro_test, h = ht_dendro)$lower[[4]],
#      main="Second branch of lower tree with cut at h=25")
# 
#  fwrite(order_su, "./outputs/working_su.csv")
```
The cut-offs representing each level are approximate. Really only the lowest level c
```{r cut cluster }
 assocs <- cutree(as.hclust(ss_clst),h =.25)
alliances <- cutree(as.hclust(ss_clst), h=.75)
 suborders <- cutree(as.hclust(ss_clst), h=.85)
 orders <- cutree(as.hclust(ss_clst), h=.99)
 k = length(unique(assocs))
 k2 <- length(unique(alliances))
 k3 <- length(unique(suborders))
 k4 <- length(unique(orders))
# ss_dend <- fviz_dend (ss_clst , k=k4, color_labels_by_k = T, lwd = .5, rect = T, cex = .5, horiz = T) %>%    plot()
  # 

 assocs_su <- as.data.frame(assocs) %>% rownames_to_column("SiteUnit")# %>% mutate(working = paste0("working-", orders))
 assocs_su <- left_join(su, assocs_su, by = "SiteUnit")
 fwrite(assocs_su, "./outputs/Assoc_w_plots_su.csv")
 groups <- assocs_su %>% column_to_rownames("PlotNumber") %>% dplyr::select(assocs) %>% t 
 
assocs_su <- as.data.frame(assocs) %>% rownames_to_column("SiteUnit") %>% mutate(assocs = paste0("workgrp-", assocs))
assocs_SS <- left_join(su, assocs_su, by = "SiteUnit") %>% dplyr::select (SiteUnit, assocs) %>% distinct %>% drop_na
# %>% arrange(desc(PlotNumber))
alliance_su <- as.data.frame(alliances) %>% rownames_to_column("SiteUnit") %>% mutate(alliances = paste0("alliance-", alliances))
suborder_su <- as.data.frame(suborders) %>% rownames_to_column("SiteUnit") %>% mutate(suborders = paste0("suborder-", suborders))
newsuborder <- unique(suborder_su$suborders)
order_su <- as.data.frame(orders) %>% rownames_to_column("SiteUnit") %>% mutate(orders = paste0("order-", orders))
#order_su <- order_su %>% mutate(orders = ifelse((orders %in% fewplots$x), "Unplaced", orders))
#neworder <- unique(order_su$orders)

Hier_su <- left_join(su, assocs_su, by = "SiteUnit") %>% 
left_join(alliance_su, by = "SiteUnit") %>% 
  left_join(suborder_su, by = "SiteUnit") %>% 
  left_join(order_su, by = "SiteUnit") %>% distinct() %>% drop_na() %>% select(-row_names)

#AllUnits_su <- left_join(su, Hier_su , by = "SiteUnit") %>% distinct
#fwrite(AllUnits_su, "./outputs/Coast_AllForestUnits_su.csv")
working_count <- length(unique(Hier_su$assocs))
ss_count <- length(unique(Hier_su$SiteUnit))
fwrite(Hier_su, "./outputs/Coast_ForestClusteredUnits_su.csv")
```

Build a Vpro hierachy style to be copied into a copy of the sample hier and then save as to a real hierarchy where the ID can be turned back into an autonumber for proper working of Vpro hierarchy functions

```{r Convert to hierarchy}
Hier.new <- Hier_su %>% mutate(formation = "-", class = "-", suball = "-", subass = "-", facies = "-", assocs = "-", working = "-") %>%   
  dplyr::select(formation, class, orders, suborders,alliances, suball, assocs, subass, facies, working, SiteUnit) %>% 
  distinct()
Hier.new  <- Hier.new  %>% mutate(class = ifelse(is.na(working), "unplaced", "placed" ))
#Hier.new <- fread("./outputs/BECv13_Hierarchy_Matrix.csv")
 levelNames <- c("formation", "class", "orders", "suborders", "alliances", "suball", "assocs", "subass", "facies", "working", "SiteUnit")
Hier.new <- as.data.table(Hier.new)
testReverse <- tableToTree(hierWide = copy(Hier.new),levelNames)
testReverse$Parent <- ifelse(testReverse$Parent == 1, "", testReverse$Parent)
fwrite(testReverse, "./outputs/SIFRForestHierarchyTreev2.csv")
```

####______________SOME OTHER TESTS_______________________

Build assocs in two steps. Build an initial very similar set and then cluster the initial cluster again
```{r cutree}

 assocs <- cutree(as.hclust(ss_clst),h =.25)
 assoc_su <- as.data.frame(assocs) %>% rownames_to_column("SiteUnit") %>% mutate(assocs = paste0("working-", assocs))
assoc_SS <- left_join(su_reduced, assoc_su, by = "SiteUnit") %>% dplyr::select (SiteUnit, assocs) %>% distinct %>% drop_na
count_assocs <- length(unique(assoc_SS$assocs))
```
Repeat cluster with first round working and look for similar

```{r rebuild the hierarchy using the new assoc as the root}
su2 <- left_join(su_reduced, assoc_su) %>% dplyr::rename(siteseries = SiteUnit, SiteUnit = assocs)

source("./_functions/_create_su_vegdata.R")
vegsum <- create_su_vegdata(vegdata, su2)

source("./_functions/_create_analysis_vegsum.R")
veg_anal <- create_analysis_vegsum (vegsum, importance = 0.5, constancy = 60, minplots = 5)

veg.dat <-
  veg_anal %>% dplyr::select(SiteUnit, Species, spp_importance) %>% matrify()
# test_jac.vegan <- vegan::vegdist(test.data, method = "sorensen", binary = F) %>% as.matrix
#test_jac <- philentropy::distance(test.data, method = "jaccard")
test_jac <- proxy::dist(veg.dat, method="ejaccard", diag=TRUE, upper = TRUE)
ss_clst2 = agnes(test_jac, diss = TRUE, stand = FALSE,
      method = "complete")
 assocs2 <- cutree(as.hclust(ss_clst2),h =.25)
 assoc_SS2 <- as.data.frame(assocs2) %>% rownames_to_column("SiteUnit") %>% mutate(assocs2 = paste0("assoc-", assocs2)) %>% dplyr::rename(assocs = SiteUnit)
 assoc_su2 <- left_join(assoc_su, assoc_SS2, by = "assocs")
 su3 <- su2 %>% dplyr::rename(assocs = SiteUnit, SiteUnit = siteseries) %>% dplyr::select(-assocs)
assoc_su2 <- left_join(su3, assoc_su2, by = "SiteUnit") %>% dplyr::select (PlotNumber, SiteUnit, assocs, assocs2) %>% dplyr::rename(working = assocs, assocs = assocs2) %>% distinct %>% drop_na
assoc_SS2 <- assoc_su2 %>% dplyr::select (-PlotNumber) %>% distinct %>% drop_na
count_assocs2 <- length(unique(assoc_SS2$assocs))
#Hier_su <- Hier_su %>% dplyr::select(-assocs)
assoc_su2 <- assoc_su2 %>% dplyr::select(SiteUnit, assocs)
Hier_su2 <- left_join(Hier_su, assoc_su2, by = "SiteUnit")

# working <- cutree(as.hclust(ss_clst),h =.25)
# alliances <- cutree(as.hclust(ss_clst), h=.5)
#  suborders <- cutree(as.hclust(ss_clst), h=.80)
#  orders <- cutree(as.hclust(ss_clst), h=1)
#   # 
# alliance_su <- as.data.frame(alliances) %>% rownames_to_column("SiteUnit") %>% mutate(alliances = paste0("alliance-", alliances))
# suborder_su <- as.data.frame(suborders) %>% rownames_to_column("SiteUnit") %>% mutate(suborders = paste0("suborder-", suborders))
# newsuborder <- unique(suborder_su$suborders)
# order_su <- as.data.frame(orders) %>% rownames_to_column("SiteUnit") %>% mutate(orders = paste0("order-", orders))
# #order_su <- order_su %>% mutate(orders = ifelse((orders %in% fewplots$x), "Unplaced", orders))
# #neworder <- unique(order_su$orders)

Hier_su2 <- left_join(assoc_su2, alliance_su, by = "SiteUnit") %>% 
  left_join(suborder_su, by = "SiteUnit") %>% 
  left_join(order_su, by = "SiteUnit") %>% distinct() %>% drop_na()

fwrite(Hier_su, "./outputs/CoastForest2stepClusteredUnits_su.csv")
```

```{r rebuild the hierarchy using the new assoc as the root}
Hier.new <- Hier_su2 %>% mutate(formation = "-", class = "-", suball = "-", subass = "-", facies = "-") %>%   
  dplyr::select(formation, class, orders, suborders,alliances, suball, assocs, subass, facies, working, SiteUnit) %>% 
  distinct()
Hier.new  <- Hier.new  %>% mutate(class = ifelse(is.na(working), "unplaced", "placed" ))
#Hier.new <- fread("./outputs/BECv13_Hierarchy_Matrix.csv")
 levelNames <- c("formation", "class", "orders", "suborders", "alliances", "suball", "assocs", "subass", "facies", "working", "SiteUnit")
Hier.new <- as.data.table(Hier.new)
testReverse <- tableToTree(hierWide = copy(Hier.new),levelNames)
testReverse$Parent <- ifelse(testReverse$Parent == 1, "", testReverse$Parent)
fwrite(testReverse, "./outputs/CoastForestHierarchyTree2.csv")


```

Build groups iteratively merging pairs of units with the minimum distance (and <.20) and add an assoc column to su table for each of the site series in the pair with an assoc number. Then recalculating veg_anal from the list of units a Set all site series to equivalent plot numbers. And iterate.
```{r novel grouping method}
# xx <- test_jac.list[ , .SD[which.min(dist)], by = iso2] %>% filter(dist < 0.2) %>% arrange(iso1)
# unit_list <- unique(xx$iso1) %>% data.frame
# xx[, ID := sequence(.N), by = iso1]
# xx[, change := shift(iso1) != iso1]
# xx[is.na(change), change := TRUE]
# xx[, res := cumsum(change)]
# 
# xx1 <- xx %>% dplyr::select(iso1, res) %>% rename(SiteUnit = 1) %>% distinct
# xx2 <- xx %>% dplyr::select(iso2, res) %>% rename(SiteUnit = 1) %>% distinct
# xx_all <- rbind(xx1, xx2)
```