tot_div_birds.Rmd

---
title: "tot_div_birds"
author: "Robert S. C. Cooke (03rcooke@gmail.com)"
date: "22/06/2020"
output: html_document
---

## Setup required packages

#### readr: read data
#### dplyr: data manipulation
#### tidyr: data manipulation
#### purrr: applying functions
#### jtools: model summary stats
#### rsq: partial r2
#### broom: model coefficients
#### arm: posterior distribution of coefficients
#### sp: spatial
#### raster: spatial
#### zoo: rolling averages
#### HDInterval: highest (posterior) density intervals

```{r setup}

if(!require("pacman")) install.packages("pacman")
pacman::p_load(readr, dplyr, tidyr, purrr, jtools, rsq, broom, arm, sp, raster, zoo, HDInterval)

```

## Raw data

# island data
# Weigelt et al., 2013 - https://www.pnas.org/content/110/38/15307

# extinct birds data
# Sayol et al., 2020

# bird distributions
# Bird species distribution maps of the world. Version 2018.1. - http://datazone.birdlife.org/species/requestdis

# grid climatic data
# WorldClim; Hijmans et al., 2005
# Elevation; Farr, 2007

# landmasses and archipelagos
# GADM version 3.6 - gadm.org/data.html

# human arrival dates
# see references in 'colz'

# extinction probabilites for possibly extinct species
# Butchart et al., 2018

## Load in preprocessed data

# Preprocessing code available upon request (03rcooke@gmail.com)

```{r}

# Data S2. Recorded prehistoric non-passerines predictor data
pred <- readr::read_csv("~/R/R_Projects/total_div/data/data_S2_pred.csv")
# Archip = Archipelago; Dist = Isolation distance; SLMP = Surrounding landmass; Elev = Elevation; Temp = Temperature; Prec = Precipitation; varT = Temperature variability; varP = Precipitation variability; SR_GAM = Archipelago plant richness; SRML = Mainland plant richness; no_isl = Number of islands; tot_area = Total area; sd_area = SD area; pa_rodents = Native rodents; hum = Human arrival; end_spp_np = Endemic non-passerines; res_eff = Research effort; rec_pre_np = Recorded prehistoric non-passerines; rec_pre_pass = Recorded prehistoric passerines; ant_spp_np = Non-passerines extant at 1500 CE (extant and historic extinct); ant_spp_pass = Passerines extant at 1500 CE (table S1)

# Data S4. Prehistoric non-passerine description dates
descr <- readr::read_csv("~/R/R_Projects/total_div/data/data_S4_descr.csv")
# species = Species name; ref.date = Reference for date; authority = Taxonomic authority; described = Year of description

# Data S5. Global passerine probability predictor data
pass <- readr::read_csv("~/R/R_Projects/total_div/data/data_S5_pass.csv")
# x = Longitude; y = Latitude; Elev = Elevation; Temp = Temperature; Prec = Precipitation; varT = Temperature variability; varP = Precipitation variability; prop_x = Surrounding landmass proportion with a neighborhood of size x; np_spp = Number of non-passerines; pass_spp = Number of passerines
  
# Data S6. Archipelago predictor data
dri_arch <- readr::read_csv("~/R/R_Projects/total_div/data/data_S6_dri_arch.csv")
# Archip = Archipelago; Elev = Elevation; Temp = Temperature; Prec = Precipitation; varT = Temperature variability; varP = Precipitation variability; prop_x = Surrounding landmass proportion with a neighborhood of size x

# Data S7. Upper bound predictor data
upp <- readr::read_csv("~/R/R_Projects/total_div/data/data_S7_upp.csv")
# x = Longitude; y = Latitude; Elev = Elevation; Temp = Temperature; Prec = Precipitation; varT = Temperature variability; varP = Precipitation variability; np_spp = Number of non-passerines; pass_spp = Number of passerines; all_spp = Number of non-passerines and passerines

# Data S8. Human arrival dates
colz <- readr::read_csv("~/R/R_Projects/total_div/data/data_S8_colz.csv")
# Archip = Region; prehist = Prehistorically settled (i.e., pre-1500 CE); cont = Continental; hum = Central human arrival date estimate (years before present); hum_low = Lower estimate of human arrival date (years before present); hum_upp = Upper estimate of human arrival date (years before present); ref = Reference; full_ref = Full reference 

# Data S9. Bird extinction dates input data
bird <- readr::read_csv("~/R/R_Projects/total_div/data/data_S9_bird.csv")
# species = Species name; common = Common name; order = Order; mod_date = 0 (prehistoric), 1 (historic), 2 (possibly extinct); insular = focal (focal archipelagos; fig. S2), insular (archipelagos), continental (continental); ext_date_min = Minimum extinction date; ext_date_max = Maximum extinction date; ext_date = Central extinction date; p_records = Extinction probability based on records for possibly extinct species; p_threats = Extinction probability based on threats for possibly extinct species; fossil_max = Maximum age of fossil; fossil_min = Minimum age of fossil; Archip = Archipelago

```

## Modelling

#### transformations and standardization

```{r}

# transformations
pred_t <- pred %>% 
  # drop unused data
  dplyr::select(-c(rec_pre_pass, ant_spp_np, ant_spp_pass)) %>% 
  dplyr::mutate_at(vars(Dist, tot_area, SR_GAM, hum, Elev), log10) %>% 
  dplyr::mutate_at(vars(sd_area, res_eff, rec_pre_np, end_spp_np), ~log10(. + 1)) %>% 
  dplyr::mutate_at(vars(no_isl), sqrt) 

# centered and scaled to zero mean and unit variance
pred_ts <- pred_t %>% 
  dplyr::mutate_at(vars(-Archip, -rec_pre_np), scale) %>% 
  dplyr::mutate_at(vars(-Archip, -rec_pre_np), as.numeric)

```

#### linear model of recorded prehistoric extinct non-passerines

```{r}

m1 <- lm(rec_pre_np ~ Dist + SLMP + Elev + Temp + varT + Prec + varP + SRML + SR_GAM + no_isl + tot_area + sd_area + pa_rodents + end_spp_np + res_eff + hum + tot_area:res_eff, data = pred_ts)

# model summary
jtools::summ(m1)

# coefficients and partial r2
coef <- rsq::rsq.partial(m1, type = "sse") %>% 
  data.frame() %>% 
  dplyr::mutate(partial.rsq_round = round(partial.rsq, digits = 2)) %>% 
  dplyr::full_join(broom::tidy(m1, conf.int = TRUE), by = c("variable" = "term"))

```

## Extrapolate values

```{r}

# research effort in new zealand
nz <- pred_ts %>% 
  dplyr::filter(Archip == "New_Zealand") %>% 
  dplyr::select(res_eff)

# 1,000 posterior draws for example (400,000 used for study)
draw <- 2000

# draw the posterior coefficients
post_coef <- arm::sim(m1, draw) %>%
  .@coef %>%
  as.data.frame()

# extrapolate for each archipelago for each posterior sample
pred_dist <- lapply(1:nrow(pred_ts), function(x) {
  
  # archipelago
  ar <- pred_ts[x,]$Archip
  
  # data for archipelago
  nd <- pred_ts %>%
    dplyr::filter(Archip == ar)
  
  # recorded prehistoric extinct non-passerines
  pre_spp_np <- nd$rec_pre_np
  
  # delta research effort (compared to New Zealand)
  delta_res <- nz - nd$res_eff
  
  # for each posterior sample
  out <- lapply(1:draw, function(b) {
    
    # total prehistoric extinct non-passerines
    tot_spp_np <- pre_spp_np + delta_res$res_eff * (post_coef[["res_eff"]][b] + post_coef[["tot_area:res_eff"]][b] * post_coef[["tot_area"]][b])
    
    # transform to natural scale
    tot_spp_np <- (10 ^ tot_spp_np) - 1
    
    # subtract recorded prehistoric extinct non-passerines
    ext_spp_np <- tot_spp_np - ((10 ^ pre_spp_np) - 1)
    
  })
  
  # extract data and tidy
  out_df <- unlist(out) %>%
    as.data.frame() %>%
    setNames(ar)
  
  # print archipelago name - progress
  print(ar)
  
  # return data
  return(out_df)
  
})

# combine data across archipelagos
pred_uni <- dplyr::bind_cols(pred_dist)

```

## Lower bound - species description curve

```{r}

# description dates for lower bound - cumulative description rate
dates_descr <- dplyr::count(descr, described) %>%
  dplyr::mutate(cumsum = cumsum(n))

# functions for modelling logisitic - Daniele Silvestro
get_mu0_logistic <- function(L, k, x0, x){
	return(L / (1 + exp(-k * (x - x0))))
}

sliding_window <- function(i, wsize = 0.5){
	new_i <- i + (runif(1) - 0.5) * wsize
	return(new_i)
}

multiplier_proposal <- function(i, d = 1.1){
	u <- runif(1)
	l <- 2 * log(d)
	m <- exp(l * (u - 0.5))
 	ii <- i * m
	hastings_ratio <- log(m)
	return(c(ii, hastings_ratio))
}

# time steps
x <- dates_descr$described - min(dates_descr$described)

# cumulative description
dat <- dates_descr$cumsum

maxDat <- max(dat)
# init parameters
epsilonA <- 2   # error (variance)
LA <- maxDat * 10 # maximum of the logistic = 'dy + L'
dyA <- 0        # shifts up or down the logistic
kA <- 2         # slope
x0A <- min(x)   # mid point

mu0 <- get_mu0_logistic(LA, kA, x0A, x) + dyA

likA <- sum(dnorm(dat, mu0, epsilonA, log = TRUE))

res <- NULL

# 6,000,000 iterations for example (2,000,000,000 iterations used in study)
# running time ~10 minutes
for (it in 0:6000000){
	L  <- LA
	k  <- kA
	x0 <- x0A
	dy <- dyA
	epsilon <- epsilonA
	hasting <- 0
	r <- runif(5)
	if (r[1] < 0.1) {
		temp = multiplier_proposal(LA, 1.2)
		L = temp[1]
		hasting = hasting + temp[2]
	}
	if (r[2] < 0.1) {
		temp = multiplier_proposal(kA)
		k = temp[1]
		hasting = hasting + temp[2]
	}
	if (r[3] < 0.1) {
		x0 = sliding_window(x0A, 5)
		if (x0 > max(x)){
			x0 = max(x)- (x0 - max(x))
		}
	}
	if (r[4] < 0.1) {
		dy = sliding_window(dyA, 2)
	}
	if (r[5] < 0.1) {
		temp = multiplier_proposal(epsilon, 1.2)
		epsilon = temp[1]
		hasting = hasting + temp[2]
	}
	mu0 = get_mu0_logistic(L, k, x0, x) + dy
	lik = sum(dnorm(dat, mu0, epsilon, log = TRUE))
	if (it %% 5000 == 0){
		res = rbind(res, c(it, likA, LA, kA, x0A, dyA, epsilonA))
	}
	if ((lik - likA + hasting) > log(runif(1))){
		LA  = L
		kA  = k
		x0A = x0
		dyA = dy
		epsilonA = epsilon
		likA = lik
	}
}

# tidy output data
res <- as.data.frame(res) %>% 
  setNames(c("it", "likA", "LA", "kA", "x0A", "dyA", "epsilonA"))

# simple plot of estimated results
plot(dates_descr$described, dat); lines(dates_descr$described, get_mu0_logistic(LA, kA, x0A, x) + dyA)

# sample from lower bound

# lower bound estimates
# posterior distribution
low_samp <- res[11:nrow(res),] %>%
  # subtract described species from describable species to get undescribed species
  dplyr::mutate(low_spp = LA - sum(dates_descr$n)) %>%
  dplyr::mutate(dist = "low")

```

## Lower bound rejection sampling

```{r}

# match sample sizes from extrapolations and lower bound
pred_uni_tot <- dplyr::sample_n(pred_uni, nrow(low_samp)) %>% 
  # total per sample across archipelagos
  dplyr::mutate(pred_tot = rowSums(.)) %>% 
  # lower bounds
  dplyr::bind_cols(dplyr::select(low_samp, low_spp)) 

# lower bound rejection sampling
pred_acc <- pred_uni_tot %>%
  # accept samples where prediction > lower bound
  dplyr::filter(pred_tot > low_spp) %>%
  # remove New Zealand
  dplyr::select(-New_Zealand)

# samples tested / samples accepted
print(nrow(pred_uni_tot)); print(nrow(pred_acc))
  
```

## Ratio non-passerines to passerines

```{r}

# transformations
pass_t <- pass %>% 
  dplyr::mutate(Elev = log10(Elev + 1))

# scales for predictions
mean_scale <- sapply(dplyr::select(pass_t, Elev, Temp, Prec, varT, varP, prop_5), mean)
sd_scale <- sapply(dplyr::select(pass_t, Elev, Temp, Prec, varT, varP, prop_5), sd)

# centered and scaled to zero mean and unit variance
pass_ts <- pass_t %>% 
  dplyr::mutate_at(vars(Elev, Temp, Prec, varT, varP, dplyr::starts_with("prop")), ~as.numeric(scale(.)))

# non-spatial model
pr1 <- glm(cbind(pass_spp, np_spp) ~ Elev + Temp + Prec + varT + varP + prop_5, family = binomial(logit), data = pass_ts)

# model summary
jtools::summ(pr1)

# AICc
qpcR::AICc(pr1)

# RAC model (autocovariate derived from residuals of model with environmental predictors)

# xy coordinates
xy <- cbind(pass_ts$x, pass_ts$y)
xyz <- cbind(xy, rep(NA, nrow(xy)))

# setup raster
rast <- raster::rasterFromXYZ(xyz, res = c(113000, 113000), crs = sp::CRS("+proj=eck4 +lon_0=0 +x_0=0 +ellps=WGS84 +units=m +no_defs"))

# xy coordinates and residuals
xy_residuals <- cbind(xy, resid(pr1))

# raster of residuals
rast[raster::cellFromXY(rast, xy_residuals)] <- xy_residuals[,3]

# first-order neighborhood
focal_pass_rast <- raster::focal(rast, matrix(1, nrow = 3, ncol = 3), fun = mean, na.rm = TRUE, pad = TRUE)

# extract mean residuals from neighborhood
focal_pass <- raster::extract(focal_pass_rast, xy)

# add residuals autocovariate to predictor data
pass_ts_spa <- cbind(pass_ts, focal_pass)

# RAC model
pr1_rac <- glm(cbind(pass_spp, np_spp) ~ Elev + Temp + Prec + varT + varP + prop_5 + focal_pass, family = binomial(logit), data = pass_ts_spa)

# model summary
jtools::summ(pr1_rac)

# AICc
qpcR::AICc(pr1_rac)

## predict for focal archipelagos

# transform and tidy
dri_arch_t <- dri_arch %>% 
  # transform elevation
  dplyr::mutate(Elev = log10(Elev + 1)) %>% 
  # add residuals autocovariate
  dplyr::mutate(focal_pass = median(focal_pass)) %>% 
  # remove additional surrounding landmass proportion data
  dplyr::select(-c(prop_3, prop_7:prop_19))

# scale data to match model data
dri_arch_ts <- scale(dplyr::select(dri_arch_t, -c(Archip, focal_pass)), center = mean_scale, scale = sd_scale) %>% 
  as.data.frame() %>% 
  # rejoin unscaled columns
  dplyr::bind_cols(dplyr::select(dri_arch_t, Archip, focal_pass), .)

# predict passerine probability
pass_pred <- dplyr::bind_cols(dri_arch_ts, data.frame(pass_prob = predict(pr1_rac, dri_arch_ts, type = "response"))) %>% 
  dplyr::group_by(Archip) %>% 
  dplyr::summarise(pass_prob_m = median(pass_prob)) 

# function for probability based rounding
prob_ro <- function(ro) {
  ifelse(floor(ro) + runif(n(), 0, 1) < ro, floor(ro) + 1, floor(ro))
}

# function to calculate the number of passerines using the number of non-passerines and the passerine probability
# minimum bernoulli trial needed to produce the observed number of 'successes' (non-passerines)
# arbitrary extreme maximum of 100,000 passerines
pass_bern <- function(pass_prob, np_spp, ...) {
  min(which(cumsum(rbinom(100000, 1, (1 - pass_prob))) == np_spp))
}

# function to estimate the number of passerines for each archipelago
pp_df <- function(x) {
  
  # get archipelago name
  a <- colnames(pred_acc)[[x]]
  
  # samples for archipelago
  samp <- dplyr::select(pred_acc, starts_with(a))
  
  # passerines
  out <- samp %>% 
    setNames(c("ext_np_spp_raw")) %>%
    # passerine probability
    dplyr::mutate(pass_prob = dplyr::filter(pass_pred, Archip == a)$pass_prob_m) %>% 
    # predicted non-passerines
    dplyr::mutate(ext_np_spp = prob_ro(ext_np_spp_raw)) %>% 
    # prehistoric non-passerines
    dplyr::mutate(rec_pre_np = dplyr::filter(pred, Archip == a)$rec_pre_np) %>% 
    dplyr::mutate(rec_pre_pass = dplyr::filter(pred, Archip == a)$rec_pre_pass) %>%
    # extant non-passerines
    dplyr::mutate(ant_spp_np = dplyr::filter(pred, Archip == a)$ant_spp_np) %>% 
    # extant passerines
    dplyr::mutate(ant_spp_pass = dplyr::filter(pred, Archip == a)$ant_spp_pass) %>% 
    # all non-passerines
    dplyr::mutate(np_spp = ext_np_spp + rec_pre_np + ant_spp_np) %>% 
    # bernoulli all species
    dplyr::mutate(all_spp = purrr::pmap_int(., pass_bern)) %>% 
    # all passerines
    dplyr::mutate(pass_spp = all_spp - np_spp) %>% 
    # undiscovered prehistoric extinct passerines
    dplyr::mutate(ext_pass_spp = pass_spp - (rec_pre_pass + ant_spp_pass)) %>% 
    dplyr::mutate(ext_pass_spp = ifelse(ext_pass_spp < 0, 0, ext_pass_spp)) %>% 
    # add archipelago name
    dplyr::mutate(Archip = a) %>% 
    # rearrange columns
    dplyr::select(Archip, everything())
  
  # print archipelago - progress
  print(a)
  
  # return data
  return(out)
  
}

# undiscovered prehistoric non-passerines and passerines per archipelago
# running time ~20 minutes
pp <- lapply(1:nrow(pass_pred), pp_df)

```

## Upper bound

```{r}

# transformations
upp_t <- upp %>% 
  dplyr::mutate(Elev = log10(Elev + 1))

# scales for predictions
mean_scale <- sapply(dplyr::select(upp_t, Elev, Temp, Prec, varT, varP), mean)
sd_scale <- sapply(dplyr::select(upp_t, Elev, Temp, Prec, varT, varP), sd)

# centered and scaled to zero mean and unit variance
upp_ts <- upp_t %>% 
  dplyr::mutate_at(vars(Elev, Temp, Prec, varT, varP), ~as.numeric(scale(.)))

# non-spatial model
um1 <- lm(all_spp ~ Elev + Temp + Prec + varT + varP, data = upp_ts)

# model summary
jtools::summ(um1)

# AICc
qpcR::AICc(um1)

# RAC model (autocovariate derived from residuals of model with environmental predictors)

# xy coordinates
xy <- cbind(upp_ts$x, upp_ts$y)
xyz <- cbind(xy, rep(NA, nrow(xy)))

# setup raster
rast <- raster::rasterFromXYZ(xyz, res = c(113000, 113000), crs = sp::CRS("+proj=eck4 +lon_0=0 +x_0=0 +ellps=WGS84 +units=m +no_defs"))

# xy coordinates and residuals
xy_residuals <- cbind(xy, resid(um1))

# raster of residuals
rast[raster::cellFromXY(rast, xy_residuals)] <- xy_residuals[,3]

# first-order neighborhood
focal_upp_rast <- raster::focal(rast, matrix(1, nrow = 3, ncol = 3), fun = mean, na.rm = TRUE, pad = TRUE)

# extract mean residuals from neighborhood
focal_upp <- raster::extract(focal_upp_rast, xy)

# add residuals autocovariate to predictor data
upp_ts_spa <- cbind(upp_ts, focal_upp)

# RAC model
um1_rac <- lm(all_spp ~ Elev + Temp + Prec + varT + varP + focal_upp, data = upp_ts_spa)

# model summary
jtools::summ(um1_rac)

# AICc
qpcR::AICc(um1_rac)

## predict for focal archipelagos

# tidy data
dri_arch_t <- dri_arch_t %>% 
  # add residuals autocovariate
  dplyr::mutate(focal_upp = median(focal_upp)) %>% 
  # remove additional surrounding landmass proportion data
  dplyr::select(-c(prop_5, focal_pass))

# scale data to match model data
dri_arch_ts <- scale(dplyr::select(dri_arch_t, -c(Archip, focal_upp)), center = mean_scale, scale = sd_scale) %>% 
  as.data.frame() %>% 
  # rejoin unscaled columns
  dplyr::bind_cols(dplyr::select(dri_arch_t, Archip, focal_upp), .)

# posterior draws - match predictions
draw <- nrow(pred_acc)

# number of focal archipelagos
narch <- nrow(dri_arch_ts)

# matrix of predictors
matp <- cbind(rep(1, narch), as.matrix(dplyr::select(dri_arch_ts, Elev:varP, focal_upp)))

# draw the posterior coefficients
post_um <- arm::sim(um1_rac, draw)

# build array to hold results
arrp <- array(NA, c(draw, narch))

# use matrix multiplication to fill array
for (s in 1:draw) {
  arrp[s,] <- rnorm(narch, matp %*% post_um@coef[s,], post_um@sigma[s])
}

# summarize median per archipelago per posterior draw
upp_samp <- as.data.frame(t(arrp)) %>%
  dplyr::bind_cols(dplyr::select(dri_arch_ts, Archip), .) %>%
  dplyr::group_by(Archip) %>%
  dplyr::summarise_if(is.numeric, median) %>%
  # wide to long
  tidyr::gather(key = "sample", value = "upp", contains("V"))

# adjust for archipelago area
upp_samp_area <- upp_samp %>%
  # join archipelago area
  dplyr::left_join(dplyr::select(pred, Archip, tot_area), by = "Archip") %>%
  # log10 transform area - cell area 12,769 km^2
  dplyr::mutate(gce = log10(tot_area/12769)) %>%
  # species-area relationship
  # log10(S) = log10(c) + zlog10(A)
  dplyr::mutate(upp_area = 10^(log10(upp) + (0.25 * (gce))))

# upper bound rejection sampling

# function for upper bound rejection sampling
upp_rej <- function(x) {

  # get archipelago name
  a <- unique(pp[[x]]$Archip)

  # match sample sizes
  up_m <- upp_samp_area %>%
    dplyr::filter(Archip == a) %>%
    dplyr::select(upp_area) %>%
    dplyr::sample_n(nrow(pp[[x]]))

  # upper bound
  out <- pp[[x]] %>%
    dplyr::bind_cols(up_m) %>%
    # accept samples where prediction < upper bound
    dplyr::filter(all_spp < upp_area)

  # print archipelago name - progress
  print(a)

  # return data
  return(out)

}

# run upper bound rejection sampling function
upp_out <- lapply(1:length(pp), upp_rej)

# percentage of accepted samples
comp <- data.frame(arch = pass_pred$Archip, input = sapply(pp, nrow), output = sapply(upp_out, nrow)) %>%
  dplyr::mutate(perc = ((input - output) / input) * 100)

# how many samples to downsample to
down <- 10

# downsample to 10 samples (1,000 used in study)
upp_ds <- lapply(1:length(upp_out), function(x) {
  out <- dplyr::sample_n(upp_out[[x]], down) %>%
    dplyr::mutate(run = 1:down)
})

# collapse into single dataframe
ex <- dplyr::bind_rows(upp_ds) %>%
  dplyr::select(Archip, run, ext_np_spp, ext_pass_spp, rec_pre_np, rec_pre_pass, ant_spp_np, ant_spp_pass)

# data frame for New Zealand
nz_run <- data.frame(Archip = rep("New_Zealand", down), run = 1:down, ext_np_spp = rep(0, down), ext_pass_spp = rep(0, down), rec_pre_np = rep(pred[pred$Archip == "New_Zealand",]$rec_pre_np[1], down), rec_pre_pass = rep(pred[pred$Archip == "New_Zealand",]$rec_pre_pass[1], down), ant_spp_np = rep(pred[pred$Archip == "New_Zealand",]$ant_spp_np[1], down), ant_spp_pass = rep(pred[pred$Archip == "New_Zealand",]$ant_spp_pass[1], down))

# add New Zealand to data
ex <- dplyr::bind_rows(ex, nz_run)

```

#### Extinction chronology

```{r}

# add proportional uncertainty to human arrival dates
colz_uncert <- colz %>% 
  # only archipelagos
  dplyr::filter(cont == 0) %>% 
  # average proportional uncertainty
  dplyr::mutate(hum_range = hum_low - hum_upp) %>% 
  dplyr::mutate(hum_uncert = hum_range / hum) %>% 
  # lower estimate
  dplyr::mutate(hum_low = ifelse(is.na(hum_low), hum + (mean(.$hum_uncert, na.rm = TRUE) * hum), hum_low)) %>% 
  # upper estimate
  dplyr::mutate(hum_upp = ifelse(is.na(hum_upp), hum - (mean(.$hum_uncert, na.rm = TRUE) * hum), hum_upp))

# rate - half-life of 100 years
rate <- -log(0.5)/100

# point to truncate exponential - 90% of extinctions within 332 years
trunc <- -log(1/10)/rate

# truncated exponential
# 1,000,000 random numbers (100,000,000 used in study)
exp_trunc <- rexp(1000000, rate = rate)
exp_trunc <- exp_trunc[exp_trunc < trunc]

# divide into chunks of length 100,000 (10,000,000 used in study)
m <- 100000
x <- seq_along(exp_trunc)
exp_trunc_split <- split(exp_trunc, ceiling(x / m))

# truncated exponential - Falklands only
exp_trunc_falkland <- exp_trunc_split[[6]]
exp_trunc_falkland <- exp_trunc_falkland[exp_trunc_falkland < colz_uncert[colz_uncert$Archip == "Falkland*",]$hum]
# truncate at 260 years for Falklands (up to 1950)
exp_trunc_falkland2 <- exp_trunc_split[[8]]
exp_trunc_falkland2 <- exp_trunc_falkland2[exp_trunc_falkland2 < colz_uncert[colz_uncert$Archip == "Falkland*",]$hum]
# truncate at 260 years for Falklands (up to 1950)

# visualize truncated exponential
#plot(hist(exp_trunc))

# human settlement of archipelagos
ex_da <- ex %>% 
  dplyr::left_join(colz_uncert) %>%
  # human first arrival from uniform distribution of dates - archipelago level
  dplyr::mutate(hum_unif = purrr::pmap(list(x = hum_upp, y = hum_low), ~ runif(1, .x, .y))) %>% 
  tidyr::unnest() %>% 
  dplyr::select(Archip:ant_spp_pass, hum_unif)

# undiscovered

# non-passerines
ex_da_un_np <- ex_da %>% 
  # create row for every species, keep count column
  tidyr::uncount(ext_np_spp, .remove = FALSE) %>%
  # extinction dates from truncated exponential - species level
  dplyr::mutate(pred_date = hum_unif - exp_trunc_split[[1]][1:n()]) %>% 
  # round down to integer, i.e., year
  dplyr::mutate(pred_date = floor(pred_date)) %>% 
  # distribution
  dplyr::mutate(dist = "pre_un") %>% 
  # order
  dplyr::mutate(ord = "np")

# passerines
ex_da_un_pass <- ex_da %>% 
  # create row for every species, keep count column
  tidyr::uncount(ext_pass_spp, .remove = FALSE) %>%
  # extinction dates from truncated exponential - species level
  dplyr::mutate(pred_date = hum_unif - exp_trunc_split[[3]][1:n()]) %>% 
  # round down to integer, i.e., year
  dplyr::mutate(pred_date = floor(pred_date)) %>% 
  # distribution
  dplyr::mutate(dist = "pre_un") %>% 
  # order
  dplyr::mutate(ord = "pass")

# recorded - with species names

# prehistoric extinct non-passerines
ext_pre_np_spid <- bird %>%
  dplyr::filter(insular == "focal" & !order == "Passeriformes") %>% 
  dplyr::group_by(Archip) %>% 
  dplyr::mutate(sp_id = 1:n())

ex_da_kno_np <- ex_da %>% 
  # create row for every species, keep count column
  tidyr::uncount(rec_pre_np, .remove = FALSE) %>%
  # extinction dates from truncated exponential - species level
  dplyr::mutate(pred_date = hum_unif - exp_trunc_split[[2]][1:n()]) %>% 
  # round down to integer, i.e., year
  dplyr::mutate(pred_date = floor(pred_date)) %>% 
  # distribution
  dplyr::mutate(dist = "pre_kno") %>% 
  # order
  dplyr::mutate(ord = "np") %>% 
  # add species names
  dplyr::group_by(Archip, run) %>% 
  dplyr::mutate(sp_id = 1:n()) %>% 
  dplyr::left_join(ext_pre_np_spid, by = c("Archip", "sp_id")) %>% 
  dplyr::select(Archip:ord, species, common)

# prehistoric extinct passerines
ext_pre_pass_spid <- bird %>%
  dplyr::filter(insular == "focal" & order == "Passeriformes" ) %>%
  dplyr::group_by(Archip) %>% 
  dplyr::mutate(sp_id = 1:n())

ex_da_kno_pass <- ex_da %>% 
  # create row for every species, keep count column
  tidyr::uncount(rec_pre_pass, .remove = FALSE) %>%
  # extinction dates from truncated exponential - species level
  dplyr::mutate(pred_date = hum_unif - exp_trunc_split[[4]][1:n()]) %>% 
  # round down to integer, i.e., year
  dplyr::mutate(pred_date = floor(pred_date)) %>% 
  # distribution
  dplyr::mutate(dist = "pre_kno") %>% 
  # order
  dplyr::mutate(ord = "pass") %>% 
  # add species names
  dplyr::group_by(Archip, run) %>% 
  dplyr::mutate(sp_id = 1:n()) %>% 
  dplyr::left_join(ext_pre_pass_spid, by = c("Archip", "sp_id")) %>% 
  dplyr::select(Archip:ord, species, common)

# combine datasets
ex_da <- dplyr::bind_rows(ex_da_un_np, ex_da_kno_np, ex_da_un_pass, ex_da_kno_pass)

# historic extinct found as fossils - 18 species

hist_ins <- bird %>% 
  # prehistoric extinct (actually historic) insular species
  dplyr::filter(mod_date == 0 & insular == "insular") %>% 
  # order
  dplyr::mutate(ord_pass = ifelse(order == "Passeriformes", "pass", "np")) %>% 
  # number of extinct non-passerines and passerines per archipelago
  dplyr::count(Archip, ord_pass) %>% 
  tidyr::spread(ord_pass, n, fill = 0) %>%
  # join colonization
  dplyr::left_join(dplyr::select(colz_uncert, Archip, hum_low, hum_upp), by = "Archip") %>% 
  # replicate dataframe 10 times (1,000 used in study)
  dplyr::slice(rep(1:n(), each = down)) %>% 
  # add run identifier
  dplyr::mutate(run = rep(1:down, n()/down)) %>% 
  # human first arrival from uniform distribution of dates
  dplyr::mutate(hum_unif = purrr::pmap(list(x = hum_upp, y = hum_low), ~ runif(1, .x, .y))) %>% 
  tidyr::unnest()

# non-passerines
hist_ins_np_spid <- bird %>% 
  # prehistoric extinct (actually historic) insular non-passerines
  dplyr::filter(mod_date == 0 & insular == "insular" & order != "Passeriformes") %>% 
  # convert archipelago names to underscores
  dplyr::group_by(Archip) %>% 
  dplyr::mutate(sp_id = 1:n())

hist_ins_np <- hist_ins %>% 
  # create row for every non-passerine species
  tidyr::uncount(np) %>% 
  # Falklands
  dplyr::mutate(pred_date = ifelse(!Archip == "Falkland_Islands", hum_unif - exp_trunc_split[[5]][1:n()], hum_unif - exp_trunc_falkland[1:n()])) %>% 
  # round down to integer, i.e., year
  dplyr::mutate(pred_date = floor(pred_date)) %>% 
  # distribution
  dplyr::mutate(dist = "hist_kno") %>% 
  # order
  dplyr::mutate(ord = "np") %>% 
  # add species names
  dplyr::group_by(Archip, run) %>% 
  dplyr::mutate(sp_id = 1:n()) %>% 
  dplyr::left_join(hist_ins_np_spid, by = c("Archip", "sp_id")) %>% 
  dplyr::select(Archip, run, hum_unif, pred_date, dist, ord, species, common)

# passerines
hist_ins_pass_spid <- bird %>% 
  # prehistoric extinct (actually historic) insular passerines
  dplyr::filter(mod_date == 0 & insular == "insular" & order == "Passeriformes") %>% 
  # convert archipelago names to underscores
  dplyr::group_by(Archip) %>% 
  dplyr::mutate(sp_id = 1:n())

hist_ins_pass <- hist_ins %>% 
  # create row for every passerine species
  tidyr::uncount(pass) %>% 
  # Falklands
  dplyr::mutate(pred_date = ifelse(!Archip == "Falkland_Islands", hum_unif - exp_trunc_split[[7]][1:n()], hum_unif - exp_trunc_falkland2[1:n()])) %>% 
  # round down to integer, i.e., year
  dplyr::mutate(pred_date = floor(pred_date)) %>% 
  # distribution
  dplyr::mutate(dist = "hist_kno") %>% 
  # order
  dplyr::mutate(ord = "pass") %>% 
  # add species names
  dplyr::group_by(Archip, run) %>% 
  dplyr::mutate(sp_id = 1:n()) %>% 
  dplyr::left_join(hist_ins_pass_spid, by = c("Archip", "sp_id")) %>% 
  dplyr::select(Archip, run, hum_unif, pred_date, dist, ord, species, common)

# combine datasets
hist_ins <- dplyr::bind_rows(hist_ins_np, hist_ins_pass)

# historic extinct birds with dates - 182 species
dates_hist_ext <- bird %>% 
  # historic
  dplyr::filter(mod_date == 1) %>% 
  # distribution
  dplyr::mutate(dist = "hist_kno") %>% 
  # order
  dplyr::mutate(ord = ifelse(order == "Passeriformes", "pass", "np"))

dates_hist_ext_long <- dates_hist_ext %>% 
  # replicate dataframe 10 times (1,000 used in study)
  dplyr::slice(rep(1:n(), each = down)) %>% 
  # add run identifier
  dplyr::mutate(run = rep(1:down, n()/down)) %>% 
  # add uncertainty from uniform distribution
  dplyr::mutate(ext_date = purrr::pmap(list(x = ext_date_min, y = ext_date_max), ~ runif(1, .x, .y))) %>% 
  tidyr::unnest() %>% 
  # convert to years before present
  dplyr::mutate(pred_date = floor((1950 - ext_date)))

# historic possibly extinct - 46 species
hist_pex <- bird %>%
  # possibly extinct
  dplyr::filter(mod_date == 2) %>% 
  # extinction date in years BP
  dplyr::mutate(pred_date = 1950 - ext_date) %>% 
  # average extinction probability
  dplyr::mutate(prob_ext = 1 - rowMeans(dplyr::select(., p_records, p_threats))) %>% 
  # distribution
  dplyr::mutate(dist = "hist_un") %>% 
  # order
  dplyr::mutate(ord = ifelse(order == "Passeriformes", "pass", "np"))

hist_pex_long <- hist_pex %>% 
  # replicate dataframe 10 times (1,000 used in study)
  dplyr::slice(rep(1:n(), each = down)) %>% 
  # add run identifier
  dplyr::mutate(run = rep(1:down, n()/down)) %>% 
  # binomial probability of extinction
  dplyr::mutate(pex = purrr::pmap(list(x = prob_ext), ~ rbinom(1, 1, .x))) %>% 
  tidyr::unnest() %>%
  dplyr::filter(pex == 1)

# continental prehistoric extinct
cont_pre <- bird %>% 
  # prehistoric continental extinct species
  dplyr::filter(mod_date == 0 & insular == "continental")

# continental rate - half-life of 1000 years
rate_cont <- -log(0.5)/1000

# point to truncate exponential - 90% of extinctions within 3322 years
trunc_cont <- -log(1/10)/rate_cont

# truncated exponential
# 1,000,000 random numbers (100,000,000 used in study)
exp_trunc_cont <- rexp(1000000, rate = rate_cont)
exp_trunc_cont <- exp_trunc_cont[exp_trunc_cont < trunc_cont]

# divide into chunks of length 100,000 (10,000,000 used in study)
m <- 100000
x <- seq_along(exp_trunc_cont)
exp_trunc_cont_split <- split(exp_trunc_cont, ceiling(x/m))

# all but palearctic and indo-malay
cont_pre_excl <- cont_pre %>% 
  # excluding palearctic and indo-malay
  dplyr::filter(!Archip %in% c("Palearctic", "Indo-Malay")) %>%
  # order
  dplyr::mutate(ord_pass = ifelse(order == "Passeriformes", "pass", "np")) %>% 
  # number of extinct non-passerines and passerines per archipelago
  dplyr::count(Archip, ord_pass) %>% 
  tidyr::spread(ord_pass, n, fill = 0) %>% 
  # colonization dates for continents
  dplyr::left_join(colz, by = "Archip") %>%
  # replicate dataframe 10 times (1,000 used in study)
  dplyr::slice(rep(1:n(), each = down)) %>% 
  # add run identifier
  dplyr::mutate(run = rep(1:down, n()/down)) %>% 
  # human first arrival from uniform distribution of dates - archipelago level
  dplyr::mutate(hum_unif = purrr::pmap(list(x = hum_upp, y = hum_low), ~ runif(1, .x, .y))) %>% 
  tidyr::unnest() 

# non-passerines
cont_pre_excl_np_spid <- cont_pre %>% 
  # excluding palearctic and indo-malay non_passerines
  dplyr::filter(!Archip %in% c("Palearctic", "Indo-Malay") & order != "Passeriformes") %>% 
  dplyr::group_by(Archip) %>% 
  dplyr::mutate(sp_id = 1:n())

cont_pre_excl_np <- cont_pre_excl %>% 
  # create row for every non-passerine species
  tidyr::uncount(np) %>% 
  # extinction dates from truncated exponential - species level
  dplyr::mutate(pred_date = hum_unif - exp_trunc_cont_split[[1]][1:n()]) %>% 
  # round down to integer, i.e., year
  dplyr::mutate(pred_date = floor(pred_date)) %>% 
  # distribution
  dplyr::mutate(dist = "pre_kno") %>% 
  # order
  dplyr::mutate(ord = "np") %>% 
  # add species names
  dplyr::group_by(Archip, run) %>% 
  dplyr::mutate(sp_id = 1:n()) %>% 
  dplyr::left_join(cont_pre_excl_np_spid, by = c("Archip", "sp_id"))

# passerines
cont_pre_excl_pass_spid <- cont_pre %>% 
  # excluding palearctic and indo-malay passerines
  dplyr::filter(!Archip %in% c("Palearctic", "Indo-Malay") & order == "Passeriformes") %>% 
  dplyr::group_by(Archip) %>% 
  dplyr::mutate(sp_id = 1:n())

cont_pre_excl_pass <- cont_pre_excl %>% 
  # create row for every passerine species
  tidyr::uncount(pass) %>% 
  # extinction dates from truncated exponential - species level
  dplyr::mutate(pred_date = hum_unif - exp_trunc_cont_split[[2]][1:n()]) %>% 
  # round down to integer, i.e., year
  dplyr::mutate(pred_date = floor(pred_date)) %>% 
  # distribution
  dplyr::mutate(dist = "pre_kno") %>% 
  # order
  dplyr::mutate(ord = "pass") %>% 
  # add species names
  dplyr::group_by(Archip, run) %>% 
  dplyr::mutate(sp_id = 1:n()) %>% 
  dplyr::left_join(cont_pre_excl_pass_spid, by = c("Archip", "sp_id"))

# palearctic and indo-malay
# extinction date between maximum age of fossil and 1500 CE from uniform distribution - species level
cont_pre_pal <- cont_pre %>% 
  # only palearctic and indo-malay
  dplyr::filter(Archip %in% c("Palearctic", "Indo-Malay")) %>% 
  # colonization dates for continents
  dplyr::left_join(colz, by = "Archip") %>%
  # maximum date of fossil or colonization date
  dplyr::mutate(early_date = purrr::pmap(list(x = fossil_max, y = hum_low), ~ min(.x, .y, na.rm = TRUE))) %>% 
  tidyr::unnest() %>% 
  # replicate dataframe 10 times (1,000 used in study)
  dplyr::slice(rep(1:n(), each = down)) %>% 
  # add run identifier
  dplyr::mutate(run = rep(1:down, n()/down)) %>% 
  # extinction date between maximum age of fossil and 1500 CE from uniform distribution - species level
  dplyr::mutate(hum_unif = purrr::pmap(list(x = hum_upp, y = early_date), ~ runif(1, .x, .y))) %>% 
  tidyr::unnest() %>%
  # round down to integer, i.e., year
  dplyr::mutate(pred_date = floor(hum_unif)) %>% 
  # distribution
  dplyr::mutate(dist = "pre_kno") %>% 
  # order
  dplyr::mutate(ord = ifelse(order == "Passeriformes", "pass", "np"))

# combine palearctic with other regions
cont_pre_comb <- dplyr::bind_rows(cont_pre_excl_np, cont_pre_excl_pass, cont_pre_pal)

# combine all data
ex_da_all <- ex_da %>% # predicted prehistoric unrecorded and recorded extinctions - ~2000 spp
  dplyr::bind_rows(dates_hist_ext_long) %>% # add historic extinct with extinction dates - 182 spp
  dplyr::bind_rows(hist_ins) %>% # add historic extinct only recorded from fossils - 18 spp
  dplyr::bind_rows(hist_pex_long) %>% # add possibly extinct species - ~25 spp (up to 46 spp)
  dplyr::bind_rows(cont_pre_comb) %>% # add prehistoric continental extinctions - 98 spp
  # select needed columns
  dplyr::select(Archip, run, pred_date, dist, ord, species, common) %>% 
  # BCE/CE
  dplyr::mutate(pred_date_ce = 1950 - pred_date)

# SEE data S10 for the ex_da_all produced in this study

```

#### Extinction rate

```{r}

# extinction rate through time for individual runs
ex_da_years_all <- ex_da_all %>%
  group_by(run) %>%
  # number of extinctions per year - n
  dplyr::count(pred_date) %>%
  # add in all years in sequence with zero extinctions
  tidyr::complete(pred_date = tidyr::full_seq(c(-69, 126000), period = 1), fill = list(n = 0)) %>%
  # order by year
  dplyr::arrange(-pred_date) %>%
  # cumulative extinctions per year
  dplyr::mutate(cumsum = cumsum(n)) %>%
  # total extinctions per run
  dplyr::group_by(run) %>%
  dplyr::mutate(tex = sum(n)) %>%
  # number alive per year - extant + all extinct - extinct at time T
  dplyr::mutate(alive = (10865 + tex) - cumsum) %>%
  # extinction rate per year
  dplyr::mutate(ex_rate = n / alive) %>%
  # moving average of 100 years
  dplyr::mutate(ma = zoo::rollmean(ex_rate, k = 100, fill = NA)) %>%
  # replace tiny negative numbers
  dplyr::mutate(ma = ifelse(ma < 0, 0, ma)) %>% 
  # BCE/CE
  dplyr::mutate(pred_date_ce = 1950 - pred_date)

# mean extinction rate through time across runs
ex_da_years_sum <- ex_da_years_all %>%
  # group by year
  group_by(pred_date) %>%
  # calculate mean extinction rate across runs per year
  dplyr::summarise(ma_mean = mean(ma)) %>% 
  # BCE/CE
  dplyr::mutate(pred_date_ce = 1950 - pred_date)

# SEE data S11 for the ex_da_years_sum produced in this study

# cumulative extinctions
ex_da_years_cum <- ex_da_years_all %>%
  # group by year
  group_by(pred_date) %>%
  # calculate mean extinction rate per year
  dplyr::summarise(cum_mean = mean(cumsum)) %>% 
  # BCE/CE
  dplyr::mutate(pred_date_ce = 1950 - pred_date)

# estimated species extinction dates
sp_date <- ex_da_all %>% 
  dplyr::group_by(species, common, Archip) %>% 
  dplyr::summarise_at(vars(pred_date), list(~median(.), ~sd(.))) %>% 
  dplyr::filter(!is.na(species)) %>% 
  # BCE/CE
  dplyr::mutate(median_ce = 1950 - median)

# SEE data S12 for the sp_date produced in this study

```

#### Totals

```{r}

## Example totals ##

# easy to separate different groups using dist and ord
# dist = pre_un (undiscovered prehistoric extinct); pre_kno (recorded prehistoric extinct); hist_kno (recorded historic extinct); hist_un (possibly historic extinct)
# ord = np (non-passerines); pass (passerines)

## total

total <- ex_da_all %>%
  dplyr::count(run)

# median
tot <- median(total$n)
# HDI
tot_hdi <- HDInterval::hdi(total$n, credMass = 0.95)

## undiscovered

# undiscovered extinct birds per run (1,000 runs used in study)
pre_undis <- ex_da_all %>%
  dplyr::filter(dist == "pre_un") %>% 
  dplyr::count(run)

# median
un <- median(pre_undis$n)
# HDI
un_hdi <- HDInterval::hdi(pre_undis$n, credMass = 0.95)

```