LF_10_sup_figures.Rmd

---
title: "LF_10_sup_figures"
author: "Caroline McKeon"
date: "10/19/2021"
output: word_document
---

## DATA NEEDED:

# f_oc_sum - dataframe of results from occurrence model object - created in LF_06b_oc_estimates.R
# f_pc_sum - dataframe of results from abundance model object - created in LF_06a_pc_estimates.R

# Data_ModelDF.rds  - Model dataframe with 624696 obs of 24 variables, created in LF_01_data_handling.R

# Data_08_anthromes_2_base_data_GeoTIFF/potveg.tif data from Ellis et al. (2010)
# Data_08_average_biome_clim_vars.rds created in LF_08_clim_biome_ecoregion.Rmd
# Data_08_average_ecoregion_clim_vars.rds created in LF_08_clim_biome_ecoregion.Rmd

# Data_08_sd_model_clim_data.rds created in LF_08_clim_biome_ecoregion.Rmd
# Data_08_average_biome_diffs.rds created in LF_08_clim_biome_ecoregion.Rmd
# Data_08_average_ecoregion_diffs.rds created in LF_08_clim_biome_ecoregion.Rmd

```{r include=FALSE}
library(tidyverse)
library(broom.mixed)
library(kableExtra)
library(ggplot2)
library(sjPlot) ## for the set_theme function
library(viridis)
library(ggpubr)
library(ggeffects)
library(gtools)
library(ggridges)
library(papeR)
library(kableExtra)
library(data.table)
library(cowplot)
library(ggdark)
library(raster)
knitr::opts_chunk$set(echo = FALSE, warning=FALSE, message = FALSE,
	fig.height = 8, fig.width = 12)
'%nin%' = Negate('%in%')
```

```{r plotting themes}
# set graphical parameters

# blank theme
set_theme(
  base = theme_classic(),
  axis.title.size = 1,
  axis.textsize = 1,
  legend.size = 1,
  legend.title.size = 1,
  geom.label.size = 1,
  plot.margin = margin(1,1,1,1, "pt")
)

## text editing
get_wraper <- function(width) {
  function(x) {
    lapply(strwrap(x, width = width, simplify = FALSE), paste, collapse="\n")}}
## colour palette
cb_pal <- c("#01665e", "#5ab4ac","#c7eae5", "#d8b365", 
            "#8c510a")
```


```{r}
# read in oc model
f_oc_sum <- readRDS("f_oc_sum.rds")
f_oc_sum$land_av <- NA
f_oc_sum$land_av[which(is.na(f_oc_sum$raunk_lf) & !is.na(f_oc_sum$landuse))] <- 1

f_oc_sum$landuse <-gsub(" vegetation", "", f_oc_sum$landuse, ignore.case = TRUE)
f_oc_sum$landuse <- as.character(f_oc_sum$landuse)
f_oc_sum$landuse[f_oc_sum$landuse == "Secondary (indeterminate age)"] <- "unknown Secondary"
f_oc_sum$landuse <- factor(f_oc_sum$landuse)

## Yvonne's ordering
f_oc_sum$landuse <- factor(f_oc_sum$landuse, 
                      levels = c("Primary forest", "unknown Secondary",  "Mature secondary", 
                                 "Intermediate secondary", "Young secondary",
                                 "Plantation forest", "Primary non-forest", "Pasture", "Cropland", "Urban"))

## create column for grey bars
f_oc_sum$cat <- as.numeric(f_oc_sum$landuse)
f_oc_sum$rep[f_oc_sum$landuse %in% c("Primary forest",  "Mature secondary", "Young secondary",
                                     "Primary non-forest", "Cropland")] <- "white"
f_oc_sum$rep[f_oc_sum$landuse %in% c("unknown Secondary","Intermediate secondary",
                                     "Plantation forest", "Pasture", "Urban")] <- "grey"
```

```{r, include=FALSE}
# read in pc model
f_pc_sum <- readRDS("f_pc_sum.rds")
f_pc_sum$landuse <- as.character(f_pc_sum$landuse)
f_pc_sum$landuse[f_pc_sum$clim == "species richness"] <- "Mature secondary"
f_pc_sum <- f_pc_sum[order(-f_pc_sum$raunk_lf, f_pc_sum$landuse),]
f_pc_sum$land_av <- 0
f_pc_sum$land_av[which(is.na(f_pc_sum$raunk_lf) & !is.na(f_pc_sum$landuse))] <- 1

f_pc_sum$landuse <-gsub(" vegetation", "", f_pc_sum$landuse, ignore.case = TRUE)
f_pc_sum$landuse <- as.character(f_pc_sum$landuse)
f_pc_sum$landuse[f_pc_sum$landuse == "Secondary (indeterminate age)"] <- "unknown Secondary"
f_pc_sum$landuse <- factor(f_pc_sum$landuse)

f_pc_sum$landuse <- factor(f_pc_sum$landuse)
## Yvonne's ordering
f_pc_sum$landuse <- factor(f_pc_sum$landuse, 
                      levels = c("Primary forest", "unknown Secondary", "Mature secondary",
                                 "Intermediate secondary", "Young secondary",
                                 "Plantation forest", "Primary non-forest", "Pasture", "Cropland", "Urban"))
f_pc_sum$cat <- as.numeric(f_pc_sum$landuse)
f_pc_sum$rep[f_pc_sum$landuse %in% c("Primary forest", "Mature secondary", "Young secondary",
                                     "Primary non-forest", "Cropland")] <- "white"
f_pc_sum$rep[f_pc_sum$landuse %in% c("unknown Secondary","Intermediate secondary",
                                     "Plantation forest", "Pasture", "Urban")] <- "grey"

```

```{r}
## read in model data
if(!exists("ModelDF")) {
  if(file.exists("Data_ModelDF.rds")) {
    try(ModelDF <- readRDS("Data_ModelDF.rds"))
  } else warning("ModelDF does not exist in this directory")
}

ModelDF$Predominant_habitat <-gsub(" vegetation", "", ModelDF$Predominant_habitat, ignore.case = TRUE)
ModelDF$Predominant_habitat <- as.character(ModelDF$Predominant_habitat)
ModelDF$Predominant_habitat[ModelDF$Predominant_habitat == "Secondary (indeterminate age)"] <- "unknown Secondary"
ModelDF$Predominant_habitat <- factor(ModelDF$Predominant_habitat)

## Yvonne's ordering
ModelDF$Predominant_habitat <- factor(ModelDF$Predominant_habitat, 
                      levels = c("Primary forest", "unknown Secondary",  "Mature secondary", 
                                 "Intermediate secondary", "Young secondary",
                                 "Plantation forest", "Primary non-forest", "Pasture", "Cropland", "Urban"))

## get percent cover data for sample sizes table
mydata <- ModelDF

## handle model dataframe to get just percent cover data, with species levels in the right format
levels(ModelDF$Best_guess_binomial) <- gsub(" ", "_", levels(ModelDF$Best_guess_binomial))
mydata <- ModelDF[ModelDF$Diversity_metric == "percent cover", ]
mydata <- mydata[mydata$Measurement !=0,]
mydata <- unique(mydata)
```

*Supplementary Figures for “Land use and climate affect plant occurrence and abundance across life forms”*

## whitakers biomes with my data (I had it somehwere else but can't find)
```{r}
#whitakers <- unique(ModelDF[,c("raunk_lf", "map", "mat", "pres_abs")])
#whitakers <- unique(whitakers[whitakers$pres_abs == 1,])

# ggplot(whitakers, aes(mat, map, colour = raunk_lf)) +
#   geom_jitter(position = position_dodge(0.02),
#              size = 2, alpha = 0.8) +
#     scale_color_viridis(discrete = TRUE) #+ facet_wrap( ~ raunk_lf)
  
```


## Map of data origins 

```{r fig.height=16, fig.width= 18}
#source("plot_enviro_pca.R")
biome_rast <- raster("Data_08_anthromes_2_base_data_GeoTIFF/potveg.tif")
biome_rast <- calc(biome_rast, fun=function(x){ x[x <= 0] <- NA; return(x)} )
  par(bg= "black", bty = "l", col="white", col.main = "white", col.lab = "white", col.axis = "white")
# plot(biome_rast, col=viridis_pal(option="C")(900))
#   points(full_clim[, 'Longitude'], full_clim[,'Latitude'], pch = 21, col = "green", cex = 2)
```

```{r fig.height = 8, fig.width = 10}
biome <- readRDS("Data_08_average_biome_clim_vars.rds")
biome$biome <- factor(biome$biome, 
                      levels = c("Tropical Evergreen Woodland", "Tropical Deciduous Woodland", 
                                 "Dense Shrubland", "Savanna", "Open Shrubland", "Deserts and Barren", 
                                 "Temperate Evergreen Woodland", "Temperate Deciduous Woodland", 
                                 "Grassland and Steppe", "Mixed Woodland", "Boreal Woodland", 
                                 "Tundra"))
  biome_ave <-ggplot(biome, aes(mat, map, colour = biome)) + 
  geom_boxplot(size = 4) +
  labs(colour = "Biome", x = "MAT (degrees C)", y = "MAP (mm)") + 
  scale_color_viridis(discrete = TRUE, option = "plasma") + 
    dark_theme_gray() + theme(legend.position="bottom") 
biome_ave
```


## Appendix S3 – Visualisation of data

```{r, fig.height= 8, fig.width= 12}
## plot model data
lu <- ggplot(ModelDF, aes(x=Predominant_habitat)) +
    geom_bar(aes(x=Predominant_habitat),  fill = "white", size=1, alpha = 1, show.legend = F) +
  geom_bar(data = ModelDF[ModelDF$Diversity_metric == "percent cover" & ModelDF$Measurement != 0,], 
               aes(x=Predominant_habitat),  size=1, alpha = 1, show.legend = F) +
       theme(axis.text.x = element_text(angle = 45, vjust = 0.5)) +
               scale_x_discrete(labels = get_wraper(10)) +
  labs(title = "A", x = "Land use") 

raunk <- ggplot(ModelDF, aes(x=raunk_lf)) +
    geom_bar(aes(x=raunk_lf),  fill = "white", size=1, alpha = 1, show.legend = F) +
  geom_bar(data = ModelDF[ModelDF$Diversity_metric == "percent cover" & ModelDF$Measurement != 0,], 
               aes(x=raunk_lf),  size=1, alpha = 1, show.legend = F) +
       theme(axis.text.x = element_text(angle = 45, vjust = 0.5)) +
               scale_x_discrete(labels = get_wraper(10)) +
  labs(title = "B", x = "Life form") 

lu2 <- ggplot(ModelDF[ModelDF$Diversity_metric == "percent cover" & ModelDF$Measurement != 0,],
              aes(x=Predominant_habitat)) +
    geom_bar(aes(x=Predominant_habitat), size=1, alpha = 1, show.legend = F) +
       theme(axis.text.x = element_text(angle = 45, vjust = 0.5)) +
               scale_x_discrete(labels = get_wraper(10),drop=FALSE) +
  labs(title = "C", x = "Land use") + ylim(0, 12500)
       
raunk2 <- ggplot(ModelDF[ModelDF$Diversity_metric == "percent cover" & ModelDF$Measurement != 0,], aes(x=raunk_lf)) +
    geom_bar(aes(x=raunk_lf),  size=1, alpha = 1, show.legend = F) +
       theme(axis.text.x = element_text(angle = 45, vjust = 0.5)) +
               scale_x_discrete(labels = get_wraper(10)) +
  labs(title = "D", x = "Life form") + ylim(0, 12500)
```

```{r, fig.height= 8, fig.width= 12}
lu <- lu + theme(axis.title.x = element_blank(), axis.title.y = element_blank(), 
                 axis.text.x = element_blank(),  axis.ticks.x = element_blank())
lu2 <- lu2 + theme(axis.title.y = element_blank()) 
raunk <- raunk + theme(axis.text.x = element_blank(), axis.title.x = element_blank(),
                         axis.ticks.x = element_blank(), axis.title.y = element_blank())
raunk2 <- raunk2 + theme(axis.title.y = element_blank()) 

resp <- ggarrange(lu, raunk, lu2, raunk2,  nrow = 2, ncol = 2, align = "hv", heights = c(2, 1)) + theme(plot.margin = margin(0.1,0.1,0.1,0.1, "cm"))
annotate_figure(resp, #top = text_grob("Relative proportions of model data", face = "bold", size = 20),
                left = text_grob("Count", size = 15, rot = 90))
```

## Data validation extended methods 

## Biome and Ecoregion visual summaries
```{r fig.height = 4, fig.width = 10}
biome <- readRDS("Data_08_average_biome_clim_vars.rds")

## make plot
biome_ave <-ggplot(biome, aes(mat, map, colour = biome)) + 
  geom_point(size = 4) +
  geom_pointrange(aes(ymin=map - map_var, ymax=map + map_var)) + 
  labs(colour = "Biome", x = "MAT (degrees C)", y = "MAP (mm)") + 
  scale_color_viridis(discrete = TRUE)

biome_var_ave <-ggplot(biome, aes(mat, mat_var, colour = biome)) + 
  geom_point(size = 4) +
  labs(colour = "Biome", x = "MAT (degrees C)", y = "MAT_var (degrees C)") + 
  scale_color_viridis(discrete = TRUE)

resp <-ggarrange(biome_ave, biome_var_ave, nrow = 1, ncol = 2, align = "hv", heights = c(1, 1), common.legend = TRUE, legend = "right") +
  theme(plot.margin = margin(0.1,0.1,0.1,0.1, "cm")) 
annotate_figure(resp, top = text_grob("Biome averages for climatic variables", face = "bold", size = 20))

```

```{r fig.height = 4, fig.width = 6}
eco <- readRDS("Data_08_average_ecoregion_clim_vars.rds")

## make plot
eco_ave <-ggplot(eco, aes(mat, map, colour = mat_var)) + 
  geom_point(size = 4) +
  geom_pointrange(aes(ymin=map - map_var, ymax=map + map_var)) + 
  labs(title = "Ecoregion averages for climatic variables", colour = "MAT_var", x = "MAT (degrees C)", y = "MAP (mm)") + 
  scale_color_viridis()


resp <-ggarrange(eco_ave, nrow = 1, ncol = 1) +
  theme(plot.margin = margin(0.1,0.1,0.1,0.1, "cm")) 
annotate_figure(resp, top = text_grob("Ecoregion averages for climatic variables", face = "bold", size = 20))

```

```{r fig.height = 6.5, fig.width = 9}
data <- readRDS("Data_08_sd_model_clim_data.rds")
biome_diff <- readRDS("Data_08_average_biome_diffs.rds")
eco_diff <- readRDS("Data_08_average_ecoregion_diffs.rds")
vars <- c("map", "mat", "map_var", "mat_var")

## plot comparison of average differences with 1 sd in model data
par(mfrow =c(3,4), mar=c(3,2,3,2))

for(i in names(data[,which(names(data) %in% vars)])){
  par(bty = "l")
  hist(data[,i], 
       breaks = 300,
       border = "grey", ylab = NULL, xlab = NULL,main = NULL, freq = F)
  lines(density(data[,i]))
  abline(v=mean(data[,i]),col="black", lwd = 2)
  abline(v=median(data[,i]),lty=2,col="black", lwd = 2)
  abline(v = mean(data[,i])+sd(data[,i]), lty = 1, col = "#A715AD", lwd = 2)
  abline(v = mean(data[,i])-sd(data[,i]), lty = 1, col = "#A715AD", lwd = 2)}

for(i in names(biome_diff[,which(names(biome_diff) %in% vars)])){
  par(bty = "l")
  hist(biome_diff[,i],
       breaks = 300, 
       border = "grey",ylab = NULL, xlab = NULL, main = NULL, freq = F)
  lines(density(biome_diff[,i]))
    abline(v=mean(biome_diff[,i]),col="black", lwd = 2)
  abline(v=median(biome_diff[,i]),lty=2,col="black", lwd = 2)
  abline(v = sd(data[,i]), lty = 1, col = "#A715AD", lwd = 2)}

for(i in names(eco_diff[,which(names(eco_diff) %in% vars)])){
  par(bty = "l")
  hist(eco_diff[,i],
       breaks = 300,
       border = "grey", ylab = NULL,main = NULL, xlab = paste(i), freq = F)
  lines(density(eco_diff[,i]))
    abline(v=mean(eco_diff[,i]),col="black", lwd = 2)
  abline(v=median(eco_diff[,i]),lty=2,col="black", lwd = 2)
  abline(v = sd(data[,i]), lty = 1, col = "#A715AD", lwd = 2)}

```


## within study locations
```{r}
w <- ModelDF[, c("Longitude", "Latitude", "SSBS", "Best_guess_binomial", "Study_name", "Region", "Predominant_habitat", "Country", "Biome",  "SS", "SSS", "SSB", "raunk_lf",  "map", "mat", "map_var", "mat_var",  "pres_abs")]
#ss <- unique(w[, c("SS", "Biome")])
ss <- unique(w[, c("SS", "Longitude", "Latitude")])
length(unique(ss$SS))
##  SS is unique in region, 75 countries (unique bar 2), 76 biomes (unique bar 3)
```

```{r}
x <- unique(ModelDF[,c("SS", "Region")])
table(x$Region)
```



```{r}
plot(biome_rast)
points(ss$Longitude, ss$Latitude)

```
```{r}
ss$island <- 1:length(ss$SS)
names(ss) <- c("SS", "x", "y", "island")

```

## get distances  between all locations within each study
```{r}
library(gdistance)
library(fasterize)
library(sf)
library(sp)
library(igraph)
library(lwgeom)
library(units)
library(geodist)
library(rgeos)
library(geosphere)
library(sfheaders)
library(rlist)

get_all_dist <- function(avi_labeled){
  ## calculate geodesic distance between all co-ordinates in dataframe 
  ## i.e. distance matrix for all points on island perimeters
  all_dist <- as.data.frame(geodist(avi_labeled[,which(names(avi_labeled) %in% c("x", "y"))],
                                  paired = FALSE,
                                  sequential = FALSE,
                                  pad = FALSE,
                                  measure = "geodesic"))
  ## relabel the distance matrix, so each distance is assocciated with points it is between
  all_dist <- cbind(avi_labeled$island,all_dist)
  names(all_dist) <- c("island", avi_labeled$island)
  ## get only lower triangle to avoid repitition
  all_dist[lower.tri(all_dist,diag=TRUE)] <- NA
  all_dist[all_dist==0]<-NA
  all_dist$island <- avi_labeled$island
  return(all_dist)}


all_dist <- list()
for(i in ss$SS){
  all_dist[[i]] <- get_all_dist(avi_labeled = ss[ss$SS == i,])
  }
```

## turn into dataframes
```{r}
for(j in 1:length(all_dist)){
  dists <- all_dist[[j]]
  inds <- as.data.frame(dists[,c(2, 1)])
  names(inds) <- c("dist", "spot")
  
  for (i in names(dists[which(names(dists) != "island")])){
  c <- as.data.frame(cbind(dists[,i], i))
  names(c) <- c("dist", "spot")
  inds <- rbind(inds, c)
  inds <- drop_na(inds)
  }
  all_dist[[j]] <- inds}

```

## get mean and maxium distances between sites within study
```{r}
c$mean <- 0
c$len <- 0
c <- c[0,]
names(c) <- c("ss", "max", "mean", "len")
for(j in names(all_dist)){
  d <- as.data.frame(cbind(j, max(as.numeric(as.character(all_dist[[j]]$dist)))/1000))
  d$mean <- mean(as.numeric(as.character(all_dist[[j]]$dist)))/1000
  d$len <- length(all_dist[[j]]$dist)
  names(d) <- c("ss", "max", "mean", "len")
  c <- rbind(c, d)
}

c$max <- as.numeric(as.character(c$max))
c$max[which(is.infinite(c$max))] <- 0
c$mean <- as.numeric(as.character(c$mean))
c$mean[which(is.infinite(c$mean))] <- 0

c$t <- as.numeric(as.character(c$mean))*c$len
sum(c$t)/sum(c$len)
```

```{r, fig.width= 5, fig.height=3}
par(mfrow = c(1,2))
c <- drop_na(c)
hist(c$max, breaks = 30, xlab = "maximum distance between sites in studies (km)")
abline(v = mean(c$max), col = "red")

hist(c$mean, breaks = 30, xlab = "mean distance between sites in studies (km)", ylab = "")
abline(v = mean(c$mean), col = "red")
```


## plot to explore
```{r}
par(mfrow = c(1,2))
far <- ss[ss$SS == "SH1_2013__Peri 1",]
x <- crop(biome_rast, extent(-90, -60,  -60, -40))
plot(biome_rast)
points(far$x, far$y)
plot(x)
points(far$x, far$y)
```


```{r}
par(mfrow = c(1,2))
far <- ss[ss$SS == "DL1_2012__Siebert 1",]
x <- crop(biome_rast, extent(20, 40,  -50, -10))
plot(biome_rast)
points(far$x, far$y)
plot(x)
points(far$x, far$y)
```
```{r}
par(mfrow = c(1,2))
far <- ss[ss$SS == "VB1_2013__Albrecht 1",]
x <- crop(biome_rast, extent(5, 17,  40, 50))
plot(biome_rast)
points(far$x, far$y)
plot(x)
points(far$x, far$y, cex = 0.8)
```


```{r}
x <- unique(ModelDF[, c("Longitude", "Latitude", "Region",  "Predominant_habitat","Country", "Biome", 
                 "SS", "SSS","SSBS", "SSB", #"raunk_lf", 
                 "map", "mat", "map_var", "mat_var")])

```

```{r}
for(i in unique(x$Region)){
  print(i)
  print(length(unique(x$SSBS[x$Region == i])))
}
```

```{r fig.width= 7, fig.height=4}
y <- unique(ModelDF[, c("Region",  "Predominant_habitat","Country", "Biome", 
                 "SS"#, "SSBS"
                 )])

ggplot(y, aes(x=Region)) +
    geom_bar(aes(x=Region),  fill = "white", size=1, alpha = 1, show.legend = F) +
  geom_bar(data = y, 
               aes(x=Region),  size=1, alpha = 1, show.legend = F) +
       theme(axis.text.x = element_text(angle = 45, vjust = 0.5)) +
               scale_x_discrete(labels = get_wraper(10)) +
  labs(title = "A", x = "Region", y = "Studies") 

ggplot(x, aes(x=Region)) +
    geom_bar(aes(x=Region),  fill = "white", size=1, alpha = 1, show.legend = F) +
  geom_bar(data = x, 
               aes(x=Region),  size=1, alpha = 1, show.legend = F) +
       theme(axis.text.x = element_text(angle = 45, vjust = 0.5)) +
               scale_x_discrete(labels = get_wraper(10)) +
  labs(title = "B", x = "Region", y = "Sites") 

ggplot(y, aes(x=Biome)) +
    geom_bar(aes(x=Biome),  fill = "white", size=1, alpha = 1, show.legend = F) +
  geom_bar(data = y, 
               aes(x=Biome),  size=1, alpha = 1, show.legend = F) +
       theme(axis.text.x = element_text(angle = 45, vjust = 0.5)) +
               scale_x_discrete(labels = get_wraper(10)) +
  labs(title = "A", x = "Biome") 
```


## prep data for ridge plots
```{r}
x <- unique(ModelDF[, c("Longitude", "Latitude", "Region",  "Predominant_habitat","Country", "Biome", 
                 "SS", "SSS","SSBS", "SSB", #"raunk_lf", 
                 "map", "mat", "map_var", "mat_var")])

mdf <- x %>% .[, which(names(.) %in% c("Predominant_habitat", "map", "mat", "map_var", "mat_var"))] 

covars <- c("map", "map_var",  "mat", "mat_var")

mdf3 <- data.frame()

for(i in covars){
  mdf$covars <-paste(i)
mdf2 <-mdf
mdf3 <- rbind(mdf3, mdf2)}

mdf <- mdf3
rm(mdf2, mdf3)

for(i in covars){
mdf$cont_value[mdf$covars == i] <- mdf[,i][mdf$covars == i]}

mdf <- mdf %>% .[, which(names(mdf) %nin% c( "map", "mat", "map_var", "mat_var"))]

```

```{r fig.width= 5, fig.height=3}
ggplot(mdf, aes(x=cont_value, y=Predominant_habitat)) + 
    labs(title = "Land use data across climates", y = "Land use type", x = "") +
  geom_density_ridges2(alpha = .3) + facet_grid(.~covars, scales='free') 


```