---
title: "Multivariate Analysis for Anthropogenic Pollen Indicators - Sustainability Special Issue" 
author: "Grant Snitker - New Mexico Consortium - Cultural Resource Sciences"
date: "17/04/2023"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

### Please reference this work using the following citation:
Please do not cite without the authors' permission. This is a paper in progress for a special issue in the journal *Sustainability*.
<br>

### Introduction
This RMarkdown demonstrates the development of anthropogenic pollen indicators for the Balearic Islands. Please note that some of the code chunks and their outputs are hidden to simplify the RMarkdown output. Please refer to the remaining headings for an outline of each analysis/procedure performed.

<br>

### 0.0 Preamble
```{r preamble, message = F, warning= F, echo = T}
# clear workspace
rm(list = ls())
# set seed
set.seed(123)
# load libraries
library(devtools)
library(vegan)
library(RColorBrewer)
library(raster)
library(dplyr)
library(stats)
library(dendextend)
library(ggplot2)
library(ggrepel)
library(indicspecies)
library(tibble)
library(foreach)
library(cluster)
library(scales)
library(forcats)
library(patchwork)
library(kableExtra)
library(adespatial)

euc.dist = function(x1, x2) sqrt(sum((x1 - x2) ^ 2))

```

<br>

### 1.0 Load pollen and vegetation datasets
```{r load, echo = T}
# Load Data
# pollen
pollen = read.csv("./data/Pollen/Pollen_data_all_taxa_with_olea_urtica_combined_with_sums.csv") # with the latest recombination of the pollen data
# Veg nitrophilous taxa only 
veg = read.csv("data/Vegetation/releves_Br_Bl_modif_transposed_updated.csv")

# reduce to common sites between veg and pollen
in.common = intersect(veg$Sample.number,pollen$Sample.number)
veg = veg[veg$Sample.number %in% in.common, ]
pollen = pollen[pollen$Sample.number %in% in.common, ]

# finish processing
#pollen
pollen$Veg.type = as.factor(pollen$Veg.type)
pollen.data  = pollen[,-c(1:4)] # remove site name, code, and veg type. These will be added again later to examine similarities between veg types. 
pollen.site.data  = pollen[,c(1:4)] # preserve the veg type, site name, and code
pollen.site.data$Veg.type  = droplevels(pollen.site.data$Veg.type) # drop levels for future labeling
# fix column names
columns = colnames(pollen.data)
columns.fixed = gsub(".t", '-t', columns, fixed = T)
columns.fixed.2 = gsub(".", ' ', columns.fixed, fixed = T)
columns.fixed.3 = gsub("major media", 'major/media', columns.fixed.2, fixed = T)
colnames(pollen.data) = columns.fixed.3 

#remove olea from pollen data
pollen.data = subset(pollen.data, select = - c(Olea))  
```
#### *1.1 Perform sqrt transformation on the pollen dataset for all future analyses*
```{r sqrt, echo = T}
pollen.data.sqrt = sqrt(pollen.data)
rownames(pollen.data.sqrt) = pollen.site.data$Sample.Veg.type#

# Create passive pollen observations (pollen sums)
pollen.data.sqrt$`Brassicaceae SUM` = 0.000001*pollen.data.sqrt$`Brassicaceae SUM`
pollen.data.sqrt$`Plantago SUM` = 0.000001*pollen.data.sqrt$`Plantago SUM`
pollen.data.sqrt$`Cerealia SUM` = 0.000001*pollen.data.sqrt$`Cerealia SUM`
pollen.data.sqrt$`Cistaceae SUM` = 0.000001*pollen.data.sqrt$`Cistaceae SUM`
pollen.data.sqrt$`Ericaceae SUM` = 0.000001*pollen.data.sqrt$`Ericaceae SUM`

```

<br>

### 2.0 Conduct Detrended Correspondence analysis to evaluate appropriate ordination technique after Lepš & Šmilauer (2003)
```{r DCA, echo = T}
pollen.dca = decorana(pollen.data.sqrt) # Lepš & Šmilauer (2003) suggest performing all transformations prior to DCA. 
pollen.dca
# axis length for DCA1 < 3 (DCA1 = 2.4), so linear ordination is preferred  

```
<br>

### 3.0 RDA and Environmental Variables
#### *3.1 Load, clean, and select appropriate environmental variables*
```{r RDA env cut, warning=FALSE, fig.align = "center",fig.width=10, fig.height=10, echo= F}
# load and select environmental data
env.data.raw = read.csv("./data/Env Variables/Env_variables_sustainability_selected.csv")
row.names(env.data.raw) = env.data.raw$Sample.number
env.data = env.data.raw[rownames(pollen.data.sqrt),]
env.data = env.data[,2:(length(env.data))]
env.data = as.data.frame(scale(env.data))

#conduct initial RDA
pollen.rda = rda(pollen.data.sqrt ~., data  = env.data)

# conduct vif
vif.cca(pollen.rda)

#Conduct cut by removing highly correlated variables (those with a variance inflation factor of approximately >20)
env.data.cut.20 = subset(env.data, select=-c(CSP_gHM,Dist_Agriculture_meters, AZ_CA, AZ_NumP, AZ_TE, AZ_AWMSI, AGRI_CA_1, AGRI_NumP_1, AGRI_TE_1, AGRI_AWMSI_1, FSN_CA_12, FSN_NumP_12, FSN_TE_12, FSN_AWMSI_12, Soil.Sealing_1000_sqmeters, Forest.Tree_1000_sqmeters, ForestTree_patches1000, ForestTreeP1000,Artificial_surfaces_sqmeters,Agriculture_sqmeters, Forests_sqmeters, Scrub_sqmeters, Open_bare_sqmeters, Water_marshes_sqmeters, Ocean_sqmeter, direct_human_pressure_0_4))

# # examine the result of cutting variables
 pollen.rda.cut.20 = rda(pollen.data.sqrt ~., data  = env.data.cut.20)
 vif.cca(pollen.rda.cut.20) # variables with vif above 20 are cut
```
#### *3.2 Run RDA without permutation selection of environmental variables*
```{r RDA no selection, warning=FALSE, fig.align = "center",fig.width=10, fig.height=10, echo = T}
set.seed(123)
# variables over 20
pollen.rda.20 = rda(pollen.data.sqrt ~., data  = env.data.cut.20)
stp = rda(pollen.data.sqrt ~ 1, data  = env.data.cut.20)
#plot(pollen.rda.20, main = "Non-redundant env variables, vif below 20")
anova(pollen.rda.20, by = "margin")
summary(pollen.rda.20)

```

#### *3.3 Plot RDA results*
```{r RDA perm plot, warning=FALSE, fig.align = "center",fig.width=10, fig.height=10, echo = F}
 # env variables and taxa

RDA.sites = scores(pollen.rda.20, display = "sites", scaling = 1) %>% as.data.frame() %>% rownames_to_column("Site_Code_Veg") 

RDA.env = scores(pollen.rda.20, display = "bp", scaling = 2) %>% as.data.frame()

RDA.env.herb = RDA.env[1,]
RDA.env.agro = RDA.env[2,]
RDA.env.tramp = RDA.env[3,]
RDA.env.landopen = RDA.env[4,]
RDA.env.fire = RDA.env[5,]
RDA.env.disturb = RDA.env[6,]
RDA.env.distfst = RDA.env[7,]
RDA.env.pasture = RDA.env[8,]

# Sites
ggplot() +
  geom_point(data = RDA.sites, aes(x = RDA1, y = RDA2), size = 2.5, shape = 1) +
  geom_text_repel(data = RDA.sites,aes(x = RDA1, y = RDA2,label = Site_Code_Veg, segment.alpha = 0.5, segment.size = 0.3), max.overlaps = 30,box.padding = unit(0.25, "lines"), size = 3.2) + 
  geom_text(data = RDA.env.herb, aes(x = RDA1, y = RDA2), label = "Herbivory pressure", nudge_x = .3, nudge_y = 0, size =4, fontface = 'bold') +
  geom_text(data = RDA.env.agro, aes(x = RDA1, y = RDA2), label = "Agropastoral use", nudge_x = 0,  nudge_y = .02, size =4, fontface = 'bold') +
  geom_text(data = RDA.env.tramp, aes(x = RDA1, y = RDA2), label = "Trampling", nudge_x =0.04, nudge_y = -0.02, size =4, fontface = 'bold') +
    geom_text(data = RDA.env.landopen, aes(x = RDA1, y = RDA2), label = "Landscape openness", nudge_x =-0.3, nudge_y = 0.015, size =4, fontface = 'bold') +
    geom_text(data = RDA.env.fire, aes(x = RDA1, y = RDA2), label = "Fire Occurence", nudge_x = 0, nudge_y = -0.03, size =4, fontface = 'bold') +
  geom_text(data = RDA.env.disturb, aes(x = RDA1, y = RDA2), label ="Distance to Urban Areas (m)", nudge_x = 0.2,  nudge_y = .022, size =4, fontface = 'bold') +
  geom_text(data = RDA.env.distfst, aes(x = RDA1, y = RDA2), label = "Distance to Forested Areas (m)", nudge_x =-0.05, nudge_y = -0.04, size =4, fontface = 'bold') +
    geom_text(data = RDA.env.pasture, aes(x = RDA1, y = RDA2), label = "Pasture (sq m)", nudge_x =-.08, nudge_y = .035, size =4, fontface = 'bold') +
  geom_vline(xintercept = c(0), color = "grey70", linetype = 2) +
  geom_hline(yintercept = c(0), color = "grey70", linetype = 2) +
  geom_segment(data = RDA.env, aes(x = 0, y = 0, xend = RDA1, yend = RDA2), alpha = 0.5, arrow = arrow(length = unit(0.2, "cm"))) +
  theme_bw() + theme(axis.title = element_text(size=13), axis.text = element_text(size = 11))   +
  xlim(-1.35, .9) +
  ylim(-.8, .7)+
  labs(x = "\nRDA1 (12.14%)", y = "RDA2 (6.49%)\n", title = "")
 
 # ggsave("./figures/RDA_Sites.svg", width = 9.1, height = 9.1, units = "in")
 # ggsave("./figures/RDA_Sites.pdf", width = 9.1, height = 9.1, units = "in")


# Species

RDAspecies = scores(pollen.rda.20, display = "species", scaling = "species") %>% as.data.frame() 
RDAspecies.matrix = as.matrix(RDAspecies)
euc.dist = function(x1, x2) sqrt(sum((x1 - x2) ^ 2))
origin = as.matrix(c(0,0))

RDAspecies$dist = foreach(i = 1:nrow(RDAspecies), .combine = c ) %do% euc.dist(RDAspecies[i,],origin[1,])
top.species = RDAspecies %>% arrange(desc(dist)) %>% slice(1:(length(RDAspecies$dist)*.50))
# add pollen sums for visualization
top.species = rbind(top.species, RDAspecies[c(57:61),])
top.species$species = row.names(top.species)

RDA.env = scores(pollen.rda.20, display = "bp", scaling = "species") %>% as.data.frame()

#create plot
ggplot() +
  geom_point(data = top.species, aes(x = RDA1, y = RDA2), size = 2, shape = 3 ) +
  geom_text_repel(data = top.species,aes(x = RDA1, y = RDA2,label = species, segment.alpha = 0.5, segment.size = 0.3), max.overlaps = 20,box.padding = unit(0.45, "lines"), size = 3.2) + 
  geom_text(data = RDA.env.herb, aes(x = RDA1, y = RDA2), label = "Herbivory pressure", nudge_x = .25, nudge_y = .015, size =4, fontface = 'bold') +
  geom_text(data = RDA.env.agro, aes(x = RDA1, y = RDA2), label = "Agropastoral use", nudge_x = 0,  nudge_y = .02, size =4, fontface = 'bold') +
  geom_text(data = RDA.env.tramp, aes(x = RDA1, y = RDA2), label = "Trampling", nudge_x =0.04, nudge_y = -0.02, size =4, fontface = 'bold') +
    geom_text(data = RDA.env.landopen, aes(x = RDA1, y = RDA2), label = "Landscape openness", nudge_x =-0.17, nudge_y = 0.03, size =4, fontface = 'bold') +
    geom_text(data = RDA.env.fire, aes(x = RDA1, y = RDA2), label = "Fire Occurence", nudge_x = 0, nudge_y = -0.03, size =4, fontface = 'bold') +
  geom_text(data = RDA.env.disturb, aes(x = RDA1, y = RDA2), label ="Distance to Urban Areas (m)", nudge_x = 0.19,  nudge_y = .03, size =4, fontface = 'bold') +
  geom_text(data = RDA.env.distfst, aes(x = RDA1, y = RDA2), label = "Distance to Forested Areas (m)", nudge_x =-0.05, nudge_y = -0.03, size =4, fontface = 'bold') +
    geom_text(data = RDA.env.pasture, aes(x = RDA1, y = RDA2), label = "Pasture (sq m)", nudge_x =-.08, nudge_y = .035, size =4, fontface = 'bold') +
  geom_vline(xintercept = c(0), color = "grey70", linetype = 2) +
  geom_hline(yintercept = c(0), color = "grey70", linetype = 2) +
  geom_segment(data = RDA.env, aes(x = 0, y = 0, xend = RDA1, yend = RDA2), alpha = 0.5, arrow = arrow(length = unit(0.2, "cm"))) +
  theme_bw() + theme(axis.title = element_text(size=13), axis.text = element_text(size = 11),  legend.position = c(.895,.25),legend.background = element_rect(fill="white",size=0.25, linetype="solid", colour ="black")) +
  xlim(-1.35, .9) +
  ylim(-.8, .7)+
  labs(x = "\nRDA1 (12.14%)", y = "RDA2 (6.49%)\n", title = "")

# ggsave("figures/RDA_Species.svg", width = 9.1, height = 9.1, units = "in")
# ggsave("figures/RDA_Species.pdf", width = 9.1, height = 9.1, units = "in")
```

<br>

### 4.0 Anthropogenic Indicators
#### *4.1 load the index values and prep dataframe*
```{r load index, warning=FALSE, fig.align = "center",fig.width=12, fig.height=12, echo = F}
anth.index = read.csv("./data/Anth_index/Env_variables_for_index_score_weighted.csv")

# prepare for pca
anth.index.data = anth.index[,c(2:9)]
anth.index.meta = anth.index[,c(1,10:12)]
rownames(anth.index.data) = anth.index.meta[,1]

```

#### *4.2 Assess species variable relationships using Pearson's r correlation coefficient*
##### *4.2.1 Species - env varibale relationships
```{r axis correlation, warning=FALSE, fig.align = "center",fig.width=12, fig.height=12, echo = F}

library(kableExtra)
cor.table = as.data.frame(round(t(cor(model.matrix(pollen.rda.20), predict(pollen.rda.20, rank=2))), digits = 3))

cor.table %>% 
    mutate_all(~cell_spec(., color = ifelse(. < -.5 | . > .5, "red"," black"))) %>%
    kable(escape = F, caption = "Species to Environmental variables correlations (cor > .5 or cor < -.5 in red)", col.names = c("Herbivory Pressure","Agropastoral Use",
                           "Trampling",
                           "Landscape Openess",
                           "Fire Occurrence", 
                           "Dist to Urban", 
                           "Dist to Forest", 
                           "Pasture Area")) %>%
    kable_styling()
```

#### *4.2 Create and plot PCA of Anth rank and species*
```{r plot anth scores pca, warning=FALSE, fig.align = "center",fig.width=12, fig.height=12, echo = F}

anth.rank.species = read.csv("./data/Anth_index/Anth_rank_species_low_mod_high_with_olea_with_sums_weighted.csv")
colnames(anth.rank.species)[2] = "Low/Absent"

# remove olea
anth.rank.species = subset(anth.rank.species, Sample.number != "Olea")

# prepare for pca
anth.rank.species.data = anth.rank.species[,c(2:4)]

rownames(anth.rank.species.data) = anth.rank.species[,1]

anth.rank.species.pca = rda(anth.rank.species.data)
summary(eigenvals(anth.rank.species.pca))
```

```{r nice plot,  warning=FALSE, fig.align = "center",fig.width=10, fig.height=10, echo=F}
library(ggrepel)
# ggplot version
PCAscores = scores(anth.rank.species.pca, display = "sites") %>% as.data.frame() %>% rownames_to_column("species") 
PCAvect = scores(anth.rank.species.pca, display = "species") %>% as.data.frame()
PCAvect.Low = PCAvect[1,] 
PCAvect.Mod = PCAvect[2,] 
PCAvect.High = PCAvect[3,] 

# plot PCA of -color
ggplot() +
  geom_point(data = PCAscores, aes(x = PC1, y = PC2), shape = 16, size = 2) +
  scale_color_manual(values =c("#1B9E77" ,"#D95F02" ,"#7570B3" ,"#E7298A", "#66A61E"), breaks = c(1,2,3,4,5)) +
  scale_shape(solid = FALSE) +
  geom_text_repel(data = PCAscores,aes(x = PC1, y = PC2,label = species, segment.alpha = 0.5, segment.size = 0.3), max.overlaps = 30,box.padding = unit(0.45, "lines"), size = 2.8) +
  geom_vline(xintercept = c(0), color = "grey70", linetype = 2) +
  geom_hline(yintercept = c(0), color = "grey70", linetype = 2) +
  geom_segment(data = PCAvect, aes(x = 0, y = 0, xend = PC1, yend = PC2), alpha = 0.5, arrow = arrow(length = unit(0.2, "cm"))) +
  geom_text(data = PCAvect.Low, aes(x = PC1, y = PC2, label = rownames(PCAvect.Low)), nudge_x = .1, nudge_y = 0.04, size =5) +
  geom_text(data = PCAvect.Mod, aes(x = PC1, y = PC2, label = rownames(PCAvect.Mod)), nudge_x = .07, nudge_y = -0.003, size =5) +
  geom_text(data = PCAvect.High, aes(x = PC1, y = PC2, label = rownames(PCAvect.High)), nudge_x = -.02, nudge_y = -0.04, size =5) +
  theme_bw() + theme(axis.title = element_text(size=13), axis.text = element_text(size = 11),  legend.position = c(.9,.115),legend.background = element_rect(fill="white",size=0.25, linetype="solid", colour ="black")) +
  labs(x = "\nPC1 (60.1%)",
       y = "\nPC2 (39.0%)", title = "Anthropogenic Pollen Indicators") 
```

#### *4.3 Evaluate group membership with k-means clustering*
```{r kmeans,  warning=FALSE, fig.align = "center",fig.width=10, fig.height=14.5,echo=F,results = FALSE}
set.seed(123)
anth.rank.kmeans = kmeans(anth.rank.species.data, 4)
summary(anth.rank.kmeans)
D = daisy(anth.rank.species.data)
anth.sil = silhouette(anth.rank.kmeans$cluster, D)
rownames(anth.sil) = rownames(anth.rank.species.data)
```

#### *4.4 Plots of group membership and indicator gradient*
```{r kmeans plots,  warning=FALSE, fig.align = "center",fig.width=10, fig.height=10,echo=F,results = FALSE}
######
# create plot
PCAscores = scores(anth.rank.species.pca, display = "sites") %>% as.data.frame() %>% rownames_to_column("species") 
PCAvect = scores(anth.rank.species.pca, display = "species") %>% as.data.frame()
PCAvect.Low = PCAvect[1,] 
PCAvect.Mod = PCAvect[2,] 
PCAvect.High = PCAvect[3,] 
PCAscores$Cluster = anth.rank.kmeans$cluster

# calculate distances
high.pt  =as.matrix(PCAscores[44,2:3])
anth.rank.species.data[44,]
mod.pt  =as.matrix(PCAscores[7,2:3])
anth.rank.species.data[7,]
low.pt  =as.matrix(PCAscores[3,2:3])
anth.rank.species.data[3,]
all.pts = as.matrix(PCAscores[,2:3])

high.dist = foreach(i = 1:nrow(all.pts), .combine = c ) %do% euc.dist(all.pts[i,],high.pt[1,])
mod.dist = foreach(i = 1:nrow(all.pts), .combine = c ) %do% euc.dist(all.pts[i,],mod.pt[1,])
low.dist = foreach(i = 1:nrow(all.pts), .combine = c ) %do% euc.dist(all.pts[i,],low.pt[1,])

out = rescale(c(low.dist, mod.dist, high.dist), to = c(0,1))

interval = length(out) / 3

all.dist = data.frame("low_dist" = out[1:interval], "mod_dist" = out[(interval+1):(interval*2)], "high_dist" = out[((interval*2)+1):(interval*3)])
all.dist$color = rescale(rowSums(all.dist))
rownames(all.dist) = rownames(anth.rank.species.data)

# create plot
ggplot() +
  geom_point(data = PCAscores, aes(x = PC1, y = PC2, color = as.factor(Cluster)), shape = 16, size = 2) +
  scale_color_manual(values =c("#D95F02","#1B9E77","#7570B3", "#808080"), breaks = c(3,4,1,2), labels = c("High","Moderate","Low", "Intermediate" )) +
  scale_shape(solid = FALSE) +
  geom_text_repel(data = PCAscores,aes(x = PC1, y = PC2,label = species, segment.alpha = 0.5, segment.size = 0.3), max.overlaps = 35,box.padding = unit(0.45, "lines"), size = 2.8) +
  geom_vline(xintercept = c(0), color = "grey70", linetype = 2) +
  geom_hline(yintercept = c(0), color = "grey70", linetype = 2) +
  geom_segment(data = PCAvect, aes(x = 0, y = 0, xend = PC1, yend = PC2), alpha = 0.5, arrow = arrow(length = unit(0.2, "cm"))) +
  geom_text(data = PCAvect.Low, aes(x = PC1, y = PC2, label = rownames(PCAvect.Low)), nudge_x = .1, nudge_y = -0.05, size =5) +
  geom_text(data = PCAvect.Mod, aes(x = PC1, y = PC2, label = rownames(PCAvect.Mod)), nudge_x = .07, nudge_y = -0.04, size =5) +
  geom_text(data = PCAvect.High, aes(x = PC1, y = PC2, label = rownames(PCAvect.High)), nudge_x = -.02, nudge_y = 0.05, size =5) +
  theme_bw() + theme(axis.title = element_text(size=13), axis.text = element_text(size = 11),  legend.position = c(.9,.89),legend.background = element_rect(fill="white",size=0.25, linetype="solid", colour ="black")) +
  xlim(-.9, 1.4) +
  ylim(-.9, 1)+
  labs(x = "\nPC1 (60.1%)",
       y = "\nPC2 (39.0%)", title = "", color = "Indicator Classes") 

# ggsave("figures/PCA_Clusters.svg", width = 9.2, height = 9.1, units = "in")
# ggsave("figures/PCA_Clusters.pdf", width = 9.2, height = 9.1, units = "in")

# plot PCA of -color
ggplot() +
  geom_point(data = PCAscores, aes(x = PC1, y = PC2, color = all.dist$color), shape = 16, size = 2) +
  scale_color_gradient(name = "Indicator Strength\n Index", high = "#E4191C", low = "grey75") +
  scale_shape(solid = FALSE) +
  geom_text_repel(data = PCAscores,aes(x = PC1, y = PC2,label = species, segment.alpha = 0.5, segment.size = 0.3), max.overlaps = 35,box.padding = unit(0.42, "lines"), size = 2.8) +
  geom_vline(xintercept = c(0), color = "grey70", linetype = 2) +
  geom_hline(yintercept = c(0), color = "grey70", linetype = 2) +
  geom_segment(data = PCAvect, aes(x = 0, y = 0, xend = PC1, yend = PC2), alpha = 0.5, arrow = arrow(length = unit(0.2, "cm"))) +
geom_text(data = PCAvect.Low, aes(x = PC1, y = PC2, label = rownames(PCAvect.Low)), nudge_x = .1, nudge_y = -0.05, size =5) +
  geom_text(data = PCAvect.Mod, aes(x = PC1, y = PC2, label = rownames(PCAvect.Mod)), nudge_x = .07, nudge_y = -0.04, size =5) +
  geom_text(data = PCAvect.High, aes(x = PC1, y = PC2, label = rownames(PCAvect.High)), nudge_x = -.02, nudge_y = 0.05, size =5) +
  theme_bw() + theme(axis.title = element_text(size=13), axis.text = element_text(size = 11),  legend.position = c(.9,.87),legend.background = element_rect(fill="white",size=0.25, linetype="solid", colour ="black")) +
  xlim(-.9, 1.4) +
  ylim(-.9, 1)+
  labs(x = "\nPC1 (60.1%)",
       y = "\nPC2 (39.0%)", title = "", color = "Indicator Classes") 

# ggsave("figures/PCA_Strength.svg", width = 9.2, height = 9.1, units = "in")
# ggsave("figures/PCA_Strength.pdf", width = 9.2, height = 9.1, units = "in")
```

```{r silhouette plot,  warning=FALSE, fig.align = "center",fig.width=14, fig.height=10,echo=F,results = FALSE}
######
# create silhouette plot of indicator strength values plot

#strength.values = 1 - all.dist[,1:3]
strength.values = all.dist
strength.values$kmeans_cluster = anth.rank.kmeans$cluster
strength.values$label = row.names(strength.values)

#rename all clusers to corresponding categories
strength.values$kmeans_cluster[strength.values$kmeans_cluster ==2] = "Intermediate"
strength.values$kmeans_cluster[strength.values$kmeans_cluster ==4] = "Moderate"
strength.values$kmeans_cluster[strength.values$kmeans_cluster ==1] = "Low/Absent"
strength.values$kmeans_cluster[strength.values$kmeans_cluster ==3] = "High"

# cluster high plot
high = subset(strength.values, kmeans_cluster == "High")
high.plot = ggplot(high, aes(x = fct_reorder(label, color), y = color )) +
  geom_col(fill = "#D95F02") +
  ylim(0,1) +
 # coord_flip() +
  labs(x = "", y = "Indicator Strength\n", title = "High") +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5, size  = 16), axis.title = element_text(size = 14), axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1, size = 11))
  
# cluster mod plot
mod = subset(strength.values, kmeans_cluster == "Moderate")
mod.plot = ggplot(mod, aes(x = fct_reorder(label, color), y = color )) +
  geom_col(fill = "#1B9E77") +
  ylim(0,1) +
  #coord_flip() +
  labs(x = "", y = "", title = "Moderate") +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5, size  = 16), axis.title = element_text(size = 14), 
        axis.text.y = element_blank(), axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1, size = 11))

# cluster low plot
low = subset(strength.values, kmeans_cluster == "Low/Absent")
low.plot = ggplot(low, aes(x = fct_reorder(label, color), y = color )) +
  geom_col(fill = "#7570B3") +
  ylim(0,1) +
 # coord_flip() +
  labs(x = "", y = "", title = "Low/Absent") +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5, size  = 16), axis.title = element_text(size = 14), 
        axis.text.y = element_blank(), axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1, size = 11))

# cluster int plot
int = subset(strength.values, kmeans_cluster == "Intermediate")
int.plot = ggplot(int, aes(x = fct_reorder(label, color), y = color )) +
  geom_col(fill = "#808080") +
  ylim(0,1) +
 # coord_flip() +
  labs(x = "", y = "", title = "Intermediate") +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5, size  = 16), axis.title = element_text(size = 14), 
        axis.text.y = element_blank(), axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1, size = 11))

(high.plot  |  mod.plot  |  low.plot  |  int.plot) + plot_layout(widths = c(1.2, 1.5, 1, 2.5))

# ggsave("figures/Silhouette.svg", width = 15, height = 9, units = "in")
# ggsave("figures/Silhouette.pdf", width = 15, height = 9, units = "in")
```

