2 Data

Predictions in buffer around observed location
Species for which the exact geographic coordinates of the population abundance records in the Living Planet Database are known. A buffer with 0.5 degree radius was drawn around these coordinates and abundance trends were calculated for suitable habitat inside these buffers only.

## -- Read data -- ####
years <- c(1992:2013)
source("code/FUN_interpolate.R")
lpd <- read.csv(params$lpdfile, na.strings = "NULL")
lpd$myid <- c(1:nrow(lpd))

## BUFFER
raw.hab <- read.csv(params$habfile)
raw.dens.suit <- read.csv(params$densfile_insuit)
raw.abun.suit <- read.csv(params$abunfile_insuit)
raw.dens.allLU <- read.csv(params$densfile_allLUC)
raw.abun.allLU <- read.csv(params$abunfile_allLUC)

# trends.hab <- trends.hab[complete.cases(trends.hab), ]
# trends.dens.suit  <- trends.dens.suit [complete.cases(trends.dens.suit ), ]
# trends.abun.suit  <- trends.abun.suit [complete.cases(trends.abun.suit ), ]
# trends.dens.allLU <- trends.dens.allLU[complete.cases(trends.dens.allLU), ]
# trends.abun.allLU <- trends.abun.allLU[complete.cases(trends.abun.allLU), ]

# First column contains pop.id from LPI
names(raw.hab)[1] <- "Pop.ID"
names(raw.dens.suit )[1] <- "Pop.ID"
names(raw.abun.suit )[1] <- "Pop.ID"
names(raw.dens.allLU)[1] <- "Pop.ID"
names(raw.abun.allLU)[1] <- "Pop.ID"

# Populations for which we have all predictions
mypops <- raw.hab[raw.hab$X1992 >0, "Pop.ID"]
mypops <- mypops[mypops %in% raw.dens.suit [raw.dens.suit $X1992 >0, "Pop.ID"]]
mypops <- mypops[mypops %in% raw.abun.suit [raw.abun.suit $X1992 >0, "Pop.ID"]]
mypops <- mypops[mypops %in% raw.dens.allLU[raw.dens.allLU$X1992 >0, "Pop.ID"]]
mypops <- mypops[mypops %in% raw.abun.allLU[raw.abun.allLU$X1992 >0, "Pop.ID"]]

myspecies <- lpd$Species.ID[lpd$ID %in% mypops]

raw.hab        <- raw.hab %>% filter(Pop.ID %in% mypops) %>% mutate(Species.ID = myspecies)
raw.dens.suit  <- raw.dens.suit %>% filter(Pop.ID %in% mypops) %>% mutate(Species.ID = myspecies)
raw.abun.suit  <- raw.abun.suit %>% filter(Pop.ID %in% mypops) %>% mutate(Species.ID = myspecies)
raw.dens.allLU <- raw.dens.allLU %>% filter(Pop.ID %in% mypops) %>% mutate(Species.ID = myspecies)
raw.abun.allLU <- raw.abun.allLU %>% filter(Pop.ID %in% mypops) %>% mutate(Species.ID = myspecies)

# Select data only
lpd.data <- lpd %>%
    # filter(Species.ID %in% myspecies) %>%
    dplyr::select(ID, Binomial, Species.ID, Common_name, myid, X1950:X2015)

lpd.data$Binomial <- as.character(lpd.data$Binomial)

lpd.data[, -c(1:5)] <- apply(lpd.data[, -c(1:5)], 2, as.numeric)
colnames(lpd.data) <- c(colnames(lpd.data[1:5]), 1950:2015)

lpd.interpolated <- lpd.data
lpd.interpolated[, -c(1:5)] <- t(apply(lpd.data[, -c(1:5)], 1, 
                                interpolate))

lpd.interpolated[which(lpd.interpolated[, "1992"] == 0), "1992"] <- NA # first year cannot be 0
lpd.interpolated <- lpd.interpolated[!is.na(lpd.interpolated[, "1992"]), ] #first year cannot be NA

lpd.interpolated <- droplevels(lpd.interpolated)

# subset species
lpd.interpolated <- lpd.interpolated[lpd.interpolated$Species.ID %in% myspecies, ]

4 Aggregating to multi-species index

Trends in geometric mean abundance across all species
With bootstrapped confidence intervals (10000 random samples of species)

# is.nan.data.frame <- function(x) {
# do.call(cbind, lapply(x, is.nan)) }

R <- 10000 # number of bootstrap samples

## Modelled trends
lambdas <- data.frame(year = plotyears)
indices <- data.frame(year = years)
std.err <- data.frame(year = plotyears)
indices.low <- data.frame(year = plotyears)
indices.upp <- data.frame(year = plotyears)
lambdas.low <- data.frame(year = plotyears)
lambdas.upp <- data.frame(year = plotyears)

# betas <- data.frame(row.names = lvls, "index" = NA, 
#                     "lower" = NA, "upper" = NA, rwnms = lvls)
# betas <- betas[, -which(names(betas) == "rwnms")]

for (i in c(1:length(all.lambdas))) {
    dat <- all.lambdas[[i]]
    # dat <- subset(dat, Species.ID %in% observedspecies)
    # dat <- subset(dat, ID %in% finalpops)
    names(dat)[5:length(dat)] <- years
    ii <- i+1
    
    # Average lambdas by species: average across populations
    dat.avg <- dat %>% group_by(Species.ID) %>% summarise_at(vars(as.character(years)), mean, na.rm=T)
    
    # Bootstrapping lambdas with 1000 replications: average across species
    results <- boot(data=dat.avg, statistic=bs, R=R)
    
    lambdas[1, ii] <- 0
    lambdas[2:14, ii] <- results$t0
    
    indices[1, ii] <- 1
    for (t in c(2:14)) {
            indices[t, ii] <- indices[(t-1), ii]*(10^(lambdas[t, ii]))
    }
    
    # get 95% CI of the mean lambda across speies
    lambdas.low[1, ii] <- 0
    lambdas.upp[1, ii] <- 0
    for (y in 1:13) {
            lambdas.low[y+1, ii] <- boot.ci(results, type="norm", index=y)$normal[2]
            lambdas.upp[y+1, ii] <- boot.ci(results, type="norm", index=y)$normal[3]
    }
    
    # Convert 95% CI to index scale
    indices.low[1, ii] <- 1
    indices.upp[1, ii] <- 1
    for (t in c(2:14)) {
           indices.low[t, ii] <- indices.low[(t-1), ii]*(10^(lambdas.low[t, ii]))
           indices.upp[t, ii] <- indices.upp[(t-1), ii]*(10^(lambdas.upp[t, ii]))
    }
    
    std.err.lambda <- c(0, apply(results$t, 2, sd))
    std.err[1, ii] <- 1
    for (t in c(2:14)) {
            std.err[t, ii] <- std.err[(t-1), ii]*(10^(std.err.lambda[t]))
    }
}

names(lambdas) <- nms
names(indices) <- nms
names(std.err) <- nms
names(indices.low) <- nms
names(indices.upp) <- nms
names(lambdas.low) <- nms
names(lambdas.upp) <- nms
dt.plot.obs <- indices.observed
dt.plot.mod <- indices.modelled

ggplot(data = dt.plot.mod, aes(x=year, y=value)) +
  
        # OBSERVED
        geom_point(data = dt.plot.obs, color = cls[1]) +
        geom_errorbar(data = dt.plot.obs, color = cls[1], width = 0.2,
                      aes(ymin = value-((lower-value)/(-1.96)), 
                          ymax = value+((upper-value)/(1.96)))) +
        geom_line(data = dt.plot.obs, lwd=0.5, color = cls[1]) +
        geom_ribbon(data = dt.plot.obs, lwd=0.5,
                    aes(x = year,
                        ymin = lower,
                        ymax = upper),
                    # ymin = index-1.96*se,
                    # ymax = index+1.96*se),
                    color = alpha(cls[1], 0),
                    fill = alpha(cls[1], 0.2)) +
        
        # MODELS
        geom_point(aes(colour=key, group=key)) +
        # geom_errorbar(aes(ymin=value-((lower-value)/(-1.96)),
        #                   ymax=value+((upper-value)/(1.96)),
        #                   width = 0.2, colour=key, group=key)) +
        # geom_errorbar(aes(ymin=index-se, ymax=index+se))
        geom_line(lwd=0.5, aes(colour=key, group=key)) +
        # geom_ribbon(lwd=0.5,
        #             aes(x = year,
        #                 ymin = lower,
        #                 ymax = upper,
        #                 fill = key,
        #                 group = key,
        #                 colour = key),
        #             color = alpha(1, 0)) +
        # ylim(c(0, 2.5)) +

        scale_colour_manual(values = cls, drop = F) +
        scale_fill_manual(values = alpha(cls, 0.2), drop = F) +
        facet_wrap(~key) +

        ggtitle(paste0("Change in ", length(finalpops)," populations of ",
                       length(observedspecies),
                       " African mammals"),
                "Geometric mean species abundance, relative to 1992") +
        # ylab(expression(GMA[t]~"/"~GMA[1992])) +
        ylab("Multi-species index") +
        # scale_x_continuous(breaks = c(1992:2005)) +
        theme_minimal() +
        theme(plot.title = element_text(hjust = 0, face = "bold"), legend.position = "top") +
        labs(colour = "", fill = "")

5 Validation metrics

RMSE$metric <- "RMSE"
BIAS$metric <- "Bias"

vm <- rbind(RMSE, BIAS)
 
# Visualize: Specify the comparisons you want
my_comparisons <- list( c(lbls[2], lbls[3]),
                        c(lbls[2], lbls[4]), 
                        c(lbls[2], lbls[5]),
                        c(lbls[2], lbls[5])) # add one to overcome plotting issue

ytmp <- as.data.frame(rbind(range(filter(vm, metric == "RMSE")$value),
                   range(filter(vm, metric == "Bias")$value)))
yrange <- (ytmp[, 2]-ytmp[, 1])

ymax <- c(max(filter(vm, metric == "RMSE")$value),
                    max(filter(vm, metric == "Bias")$value))
ymin <- c(min(filter(vm, metric == "RMSE")$value),
                    min(filter(vm, metric == "Bias")$value))

newymax <- round(ymax+(0.01*yrange), 1)
newymin <- round(ymin-(0.1*yrange), 1)
laby <- ymin-(0.05*yrange)

dat_text <- data.frame(metric = c(rep("RMSE", 4), rep("Bias", 4)), 
                       model = rep(levels(vm$model), 2),
                       label = c(wilcox.test(filter(RMSE, model == lbls[2])$value)$p.value,
                                 wilcox.test(filter(RMSE, model == lbls[3])$value)$p.value,
                                 wilcox.test(filter(RMSE, model == lbls[4])$value)$p.value,
                                 wilcox.test(filter(RMSE, model == lbls[5])$value)$p.value,
                                 wilcox.test(filter(BIAS, model == lbls[2])$value)$p.value,
                                 wilcox.test(filter(BIAS, model == lbls[3])$value)$p.value,
                                 wilcox.test(filter(BIAS, model == lbls[4])$value)$p.value,
                                 wilcox.test(filter(BIAS, model == lbls[5])$value)$p.value),
                       value = c(rbind(laby, laby, laby, laby)))

ranges <- data.frame(ymax = newymax,
                     # ymin = c(0.5, 0, -.75, 0.1),
                     ymin = newymin,
                     metric = c("RMSE", "Bias"),
                     model = "ESH")
ranges <- gather(ranges, "key", "value", -metric, -model) %>% dplyr::select(., -key)

tip <- yrange*0.1

Bias and RMSE

Sign of Linear Trends

linear.trends <- as.data.frame(cbind(lmObs, lmHab, 
                                     lmDens.a, lmDens.s, lmAbun.s))
# row.names(linear.trends) <- observedspecies
row.names(linear.trends) <- finalpops
colnames(linear.trends)  <- nms[2:6]
linear.pvals <- as.data.frame(cbind(lmObs.p, lmHab.p,
                                    lmDens.p.a, lmDens.p.s, lmAbun.p.s))
# row.names(linear.pvals) <- observedspecies
row.names(linear.pvals) <- finalpops
colnames(linear.pvals)  <- nms[2:6]
 
slopes <- linear.trends
# linear.trends[linear.pvals > 0.05] <- 0

signs <- sign(linear.trends)

S <- data.frame(row.names = c("Correct", 
                              "Underestimate", "Missed +", "False -", "Opposing -", 
                              "Overestimate",  "Missed -", "False +", "Opposing +"))
S[[nms[3]]] <- NA
S[[nms[4]]] <- NA
S[[nms[5]]] <- NA
S[[nms[6]]] <- NA

for (m in nms[3:6]) {
  S["Correct", m] <- sum(signs[[nms[2]]] == signs[[m]])/length(observedpops)
  S["Underestimate", m] <- sum(signs[[nms[2]]] > signs[[m]])/length(observedpops)
  S["Missed +", m] <- sum(signs[[nms[2]]] == 1 & signs[[m]] == 0)/length(observedpops)
  S["False -", m] <- sum(signs[[nms[2]]] == 0 & signs[[m]] == -1)/length(observedpops)
  S["Opposing -", m] <- sum(signs[[nms[2]]] == 1 & signs[[m]] == -1)/length(observedpops)
  S["Overestimate", m] <- sum(signs[[nms[2]]] < signs[[m]])/length(observedpops)
  S["Missed -", m] <- sum(signs[[nms[2]]] == -1 & signs[[m]] == 0)/length(observedpops)
  S["False +", m] <- sum(signs[[nms[2]]] == 0 & signs[[m]] == 1)/length(observedpops)
  S["Opposing +", m] <- sum(signs[[nms[2]]] == -1 & signs[[m]] == 1)/length(observedpops)
}

names(S) <- lbls[2:5]
 
kable(round(S[c(1,2,6), ],2), caption = "Simple classification of sign of population-level linear trends (proportion)")
Table 5.1: Simple classification of sign of population-level linear trends (proportion)
ESH Density Density in ESH Abundance in ESH
Correct 0.37 0.47 0.47 0.43
Underestimate 0.33 0.13 0.13 0.13
Overestimate 0.30 0.40 0.40 0.43

Trends significant at p<0.05 only

Table 5.2: Simple classification of sign of population-level linear trends (proportion)
ESH Density Density in ESH Abundance in ESH
Correct 0.37 0.2 0.17 0.1
Underestimate 0.33 0.3 0.33 0.4
Overestimate 0.30 0.5 0.50 0.5
Table 5.2: Full classification of sign of population-level linear trends (proportion)
ESH Density Density in ESH Abundance in ESH
Correct 0.37 0.2 0.17 0.10
Underestimate 0.33 0.3 0.33 0.40
Missed + 0.03 0.3 0.30 0.30
False - 0.07 0.0 0.00 0.07
Opposing - 0.23 0.0 0.03 0.03
Overestimate 0.30 0.5 0.50 0.50
Missed - 0.00 0.5 0.50 0.50
False + 0.03 0.0 0.00 0.00
Opposing + 0.27 0.0 0.00 0.00