Plotting population-level prediction errors (RMSE-values) from three models (habitat, density and abundance/combined) against potential explanatory variables in the LPI and mammal traits databases.

library(dplyr)
library(ggplot2)
library(gridExtra)
library(ggpubr)
library(ggsci)
library(tidyr)
cls3 <- pal_locuszoom()(4)[c(3,2,4)]
cls1 <- pal_jco()(3)[c(3)]
cls5 <- c(cls3, "#ef6548", "#1f78b4")
cls <- c(cls1, cls3, "#ef6548", "#1f78b4")
rm(cls1, cls3, cls5)

# # Load Goodness of Fit metrics
# gof <- read.csv("outputs/Africa/Africa_RMSE_vals_w_covars.csv")

# gof <- read.csv("outputs/Africa/Africa-ESACCI-SpeciesSpecific_goodness_of_fit_metrics.csv")
gof <- read.csv("outputs/Africa/Africa-ESACCI-CrossSpecies_goodness_of_fit_metrics.csv")

# Load trait data
Traits <- read.csv("data/density/fromLuca/M_Traits.csv")

# Load LPD
lpd <- read.csv("data/biodiversity/Africa/LPI_pops_Africa_spec_loc_observed.csv")

# Select data only
gof <- gof %>%
        # select_at(c("Species", "Taxid", "Popid", "RMSE_habitat", "RMSE_density", "RMSE_abundance")) %>%
        left_join(., Traits[,c('tax_id', 'Order', 'Family', 'BM_g', 'Diet')],
                  by = c('Species.ID' ='tax_id')) %>%
        left_join(., lpd,  by = c('popid' = 'ID'))  %>%
        droplevels()


## Test whether random effects are used
# lpd <- read.csv("data/biodiversity/Africa/LPI_pops_Africa_spec_loc_observed_final97pops_interpolated.csv")
# lpd$spid <- ifelse(lpd$Species.ID %in% spid$tax_id, "T", "F")
# 
# # Load id's of populations with random effects
# spid <- read.csv("data/density/fromLuca/Mammals_spid_taxonomy.csv", sep=";")
# 
# gof$spid <- ifelse(gof$Taxid %in% spid$tax_id, "T", "F")

gof$Units[gof$Units == "Groups per Km²"] <- "Number of groups per Km²"
gof$Units[gof$Units == "Mean numbers of individuals"] <- "Number of individuals"
gof$Units[gof$Units == "Individuals"] <- "Number of individuals"
gof$Units[gof$Units == "Individual"] <- "Number of individuals"
gof <- droplevels(gof) 

# gof_wide <- spread(gof, key = model, value = value)
# RMSE_wide <- filter(gof_wide, metric == "RMSE")
# RMSE_wide <- droplevels(RMSE_wide)

lbls <- c("Observed trend", "Habitat model", 
          "Density model (in ESH)", "Abundance model (in ESH)", 
          "Density model (all LU classes)", "Abundance model (all LU classes)")

gof <- filter(gof, model %in% lbls[c(2, 3, 6)])
gof <- droplevels(gof)
give.n <- function(d){
  return(data.frame(y = 1, label = paste("n =", length(d))))
}

boxplotRMSE <- function(var) { 
        h <- list()
        # RMSE from habitat models
        for (m in unique(gof$model)) {
                h[[m]] <- ggplot(data = filter(gof, metric == "RMSE" & model == m), 
                       aes(x = get(var), 
                           color=model,
                           y=value)) +
                        theme_bw() +
                        geom_hline(yintercept = mean(filter(gof, metric == "RMSE" & model == m)$value), linetype = 2) +
                        geom_boxplot() +
                        scale_color_manual(values = cls[which(lbls==m)]) +
                        stat_summary(fun.ymax = max,
                                     geom = "text", aes(y = .8, label = ..ymax..)) +
                        theme(axis.text.x = element_text(angle = 22.5, hjust = 1)) +
                        stat_compare_means(label = "p.format", method = "wilcox",cex = 3,
                                           ref.group = ".all.", 
                                           # comparisons = list(c()),
                                           label.y = 0.9) + 
                        stat_summary(fun.data = give.n, geom = "text") +
                        scale_x_discrete(name = var) +
                        ylab("RMSE") +
                        theme(legend.position = "n")
                        # facet_grid(~model)
        }
        # h
        grid.arrange(grobs = h, ncol = length(unique(gof$model)),
                     widths = rep(1, length(unique(gof$model))))
}

p-values and model fit for simple linear regression (dashed line):

## Numerical variables: body mass

lm_bm <- list()
# run models
for (m in unique(gof$model)) {
        lm_bm[[as.character(m)]] <- lm(value~BM_g, filter(gof, metric == "RMSE" & model == m))
}

h <- list()
        # RMSE from habitat models
for (m in unique(gof$model)) {
        h[[m]] <- ggplot(data = filter(gof, metric == "RMSE" & model == m), 
                         aes(x = BM_g, 
                             color=model,
                             y=value)) +
                theme_bw() +
                geom_point() +
                # geom_hline(yintercept = mean(gof$RMSE_habitat), linetype = 2) +
                scale_x_log10() +
                scale_color_manual(values = cls[which(lbls==m)]) +
                geom_smooth(method = "lm", lty = 2, se=F, col = 1, lwd = .5) +
                annotate("text", label = paste("italic(p) == " ,round(summary(lm_bm[[m]])$coefficients[2,4], 2)),
                         x = 5000, y = 0.325, hjust = 0, parse=T, vjust=0) +
                annotate("text", label = paste("R^2 == " ,round(summary(lm_bm[[m]])$r.squared, 3)),
                         x = 5000, y = 0.3, hjust = 0, parse=T, vjust=0) +
                ylab("RMSE") +
                theme(legend.position = "n")
}
        # h
        grid.arrange(grobs = h, ncol = length(unique(gof$model)), 
                     widths = rep(1, length(unique(gof$model))))

# create list of categorical variables
# see names(gof) for more options
myvars <- list("Diet", "T_biome", 
               "Red_list_category", "Criteria", "Primary_threat",
               "Data_type", "Units", "Data_source_type",  
               "Reason_for_data_collection", "Percentage")
 
# apply boxplot function to list of categorical variables
lapply(myvars, boxplotRMSE) 

## [[1]]
## TableGrob (1 x 3) "arrange": 3 grobs
##                                  z     cells    name           grob
## Habitat model                    1 (1-1,1-1) arrange gtable[layout]
## Density model (in ESH)           2 (1-1,2-2) arrange gtable[layout]
## Abundance model (all LU classes) 3 (1-1,3-3) arrange gtable[layout]
## 
## [[2]]
## TableGrob (1 x 3) "arrange": 3 grobs
##                                  z     cells    name           grob
## Habitat model                    1 (1-1,1-1) arrange gtable[layout]
## Density model (in ESH)           2 (1-1,2-2) arrange gtable[layout]
## Abundance model (all LU classes) 3 (1-1,3-3) arrange gtable[layout]
## 
## [[3]]
## TableGrob (1 x 3) "arrange": 3 grobs
##                                  z     cells    name           grob
## Habitat model                    1 (1-1,1-1) arrange gtable[layout]
## Density model (in ESH)           2 (1-1,2-2) arrange gtable[layout]
## Abundance model (all LU classes) 3 (1-1,3-3) arrange gtable[layout]
## 
## [[4]]
## TableGrob (1 x 3) "arrange": 3 grobs
##                                  z     cells    name           grob
## Habitat model                    1 (1-1,1-1) arrange gtable[layout]
## Density model (in ESH)           2 (1-1,2-2) arrange gtable[layout]
## Abundance model (all LU classes) 3 (1-1,3-3) arrange gtable[layout]
## 
## [[5]]
## TableGrob (1 x 3) "arrange": 3 grobs
##                                  z     cells    name           grob
## Habitat model                    1 (1-1,1-1) arrange gtable[layout]
## Density model (in ESH)           2 (1-1,2-2) arrange gtable[layout]
## Abundance model (all LU classes) 3 (1-1,3-3) arrange gtable[layout]
## 
## [[6]]
## TableGrob (1 x 3) "arrange": 3 grobs
##                                  z     cells    name           grob
## Habitat model                    1 (1-1,1-1) arrange gtable[layout]
## Density model (in ESH)           2 (1-1,2-2) arrange gtable[layout]
## Abundance model (all LU classes) 3 (1-1,3-3) arrange gtable[layout]
## 
## [[7]]
## TableGrob (1 x 3) "arrange": 3 grobs
##                                  z     cells    name           grob
## Habitat model                    1 (1-1,1-1) arrange gtable[layout]
## Density model (in ESH)           2 (1-1,2-2) arrange gtable[layout]
## Abundance model (all LU classes) 3 (1-1,3-3) arrange gtable[layout]
## 
## [[8]]
## TableGrob (1 x 3) "arrange": 3 grobs
##                                  z     cells    name           grob
## Habitat model                    1 (1-1,1-1) arrange gtable[layout]
## Density model (in ESH)           2 (1-1,2-2) arrange gtable[layout]
## Abundance model (all LU classes) 3 (1-1,3-3) arrange gtable[layout]
## 
## [[9]]
## TableGrob (1 x 3) "arrange": 3 grobs
##                                  z     cells    name           grob
## Habitat model                    1 (1-1,1-1) arrange gtable[layout]
## Density model (in ESH)           2 (1-1,2-2) arrange gtable[layout]
## Abundance model (all LU classes) 3 (1-1,3-3) arrange gtable[layout]
## 
## [[10]]
## TableGrob (1 x 3) "arrange": 3 grobs
##                                  z     cells    name           grob
## Habitat model                    1 (1-1,1-1) arrange gtable[layout]
## Density model (in ESH)           2 (1-1,2-2) arrange gtable[layout]
## Abundance model (all LU classes) 3 (1-1,3-3) arrange gtable[layout]