The state of equilibrium diversity was inferred from the correlation γ of diversity with speciation and extinction rates by predicting at which diversity level both rates equilibrate. Subsequently, we compared this level with the empirical diversity trajectory and calculated the 95% HPD.

# Load required packages
library(coda)
BDS <- 1:2 # Number of BDS analyses
RanSamp <- 2 # Random samples of speciation and extinction times
Fix <- paste(paste0("Fix_", BDS), rep(1:RanSamp, each = length(BDS)), sep = "_")
Fix <- sort(Fix)
N <- RanSamp * length(BDS)
M <- 90 # MCMC iterations

SameSr <- 50:350 # Same richness for all replicates

# Relationship of lambda and mu with diversity:
###############################################
DepMu <- DepLambda <- matrix(NA_real_, 
                             nrow = M * N, 
                             ncol = length(SameSr))
# Equilibrium level:
####################
EqDiv <- rep(NA_real_, M * N)

for(i in 1:N){
  # Get non-endemic diversity
  Div <- read.table("Results/Covariates/NonEnd.txt", 
                    sep = "\t", header = TRUE)
  # Get endemic diversity
  FixTmp <- read.table(paste0("Results/BD/", Fix[i], ".txt"), 
                       sep = "\t", header = TRUE)
  # Combined diversity
  FirstApp <- sort(FixTmp$ts, decreasing = TRUE)
  Increase <- rep(1, length(FirstApp))
  LastApp <- sort(FixTmp$te, decreasing = TRUE)
  Decrease <- rep(-1, length(LastApp))
  Decrease[LastApp == 0] <- 0
  Att <- data.frame(Time = c(FirstApp, LastApp), 
                    Change = c(Increase, Decrease), 
                    SR = NA_integer_)
  Att <- Att[order(Att$Time, decreasing = TRUE), ]
  Att$SR <- cumsum(Att$Change)
  Att <- Att[min(which(Att$Time <= 13.66)):which(Att$Time == 0)[1], ]
  End <- approx(Att[, c(1, 3)], xout = Div$time, method ="constant", rule = 2)$y
  Div$Div <- Div$NonEnd + End
  # RichThroughTime[i, ] <- Div$Div
  # SrThroughTime <- scale(Div$Div)
  SrScale <- (SameSr - mean(Div$Div)) / sd(Div$Div)
  ContLog  <- read.table(paste0("Results/Covariates/NoShift/", Fix[i], 
                                "_DD_0_s_13.66_13.66linSp_linEx_HP.log"), 
                                  header = TRUE, sep = "\t")
  ContLog <- ContLog[ContLog$beta == 1, ]
  ContLog <- ContLog[-c(1:round(nrow(ContLog) * 0.1)), ] # Remove 10% burnin
  DepLambdaTmp <- apply( ContLog[, c(10, 14)], 1, function(x) x[1] + x[1]*x[2]*SrScale ) * 10
  DepMuTmp <- apply( ContLog[, c(12, 16)], 1, function(x) x[1] + x[1]*x[2]*SrScale ) * 10
  DepLambdaTmp[DepLambdaTmp < 0] <- 0
  DepMuTmp[DepMuTmp < 0] <- 0
  DepLambda[ (1+(M*(i-1))):(M*i), ] <- t(DepLambdaTmp)
  DepMu[ (1+(M*(i-1))):(M*i), ] <- t(DepMuTmp)
  # DepLaTTTmp <- t(apply( ContLog[, c(10, 14)], 1, function(x) x[1] + x[1]*x[2]*SrThroughTime ) * 10)
  # DepMuTTTmp <- t(apply( ContLog[, c(12, 16)], 1, function(x) x[1] + x[1]*x[2]*SrThroughTime ) * 10)
  # Index of equilibrium diversity to look up in SameSr
  SrDepTmp <- t(DepLambdaTmp) - t(DepMuTmp)
  EqDivIndex <- sapply( 1:nrow(SrDepTmp), function(x) 
    if(any(SrDepTmp[x, ] <= 0)){min(which(SrDepTmp[x,] <= 0))} else {max(SameSr)})
  EqDiv[(1+(M*(i-1))):(M*i)] <- SameSr[EqDivIndex]
}
DepLambdaHPD <- apply(DepLambda, 2, function(x) HPDinterval(as.mcmc(x)))
DepMuHPD <- apply(DepMu, 2, function(x) HPDinterval(as.mcmc(x)))
EqDivHPD <- HPDinterval(as.mcmc(EqDiv))

par(las = 1, mar = c(5.1, 4.5, 0.1, 0.8))
plot(1, 1, xlim = c(50, 350), ylim = c(0, 6), type = "n", xaxs = "i",
     ylab = expression(paste("Rate (events ", lineage^{-1}, " ", Myr^{-1}, ")")),
     xlab = "Diversity")
rect(xleft = EqDivHPD[1], ybottom = par("usr")[3], 
     xright = EqDivHPD[2], ytop = par("usr")[4],
     col = adjustcolor("grey", alpha = 0.5), 
     border = NA)
polygon(c(SameSr, rev(SameSr)), 
        c(DepLambdaHPD[1, ], rev(DepLambdaHPD[2, ])), 
        col = adjustcolor("dodgerblue", alpha = 0.5), border = NA)
polygon(c(SameSr, rev(SameSr)), 
        c(DepMuHPD[1, ], rev(DepMuHPD[2, ])), 
        col = adjustcolor("red", alpha = 0.5), border = NA)
legend("top",  pch = 15, pt.cex = 1.5, bty = "n", 
       legend = c("Speciation", "Extinction", "Equilibrium diversity"),
       col = adjustcolor(c("dodgerblue", "red", "grey"), alpha = 0.5))
box()

Equilibrium diversity of Lake Ohrid diatoms Speciation and extinction rates depending on diatom richness and equilibrating beyond 295 species

Plattform and resources

Sys.time()
## [1] "2020-01-07 17:46:50 CET"
sessionInfo()
## R version 3.6.2 (2019-12-12)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.6 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/libblas/libblas.so.3.6.0
## LAPACK: /usr/lib/lapack/liblapack.so.3.6.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=de_DE.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=de_DE.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=de_DE.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=de_DE.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] coda_0.19-3 knitr_1.25 
## 
## loaded via a namespace (and not attached):
##  [1] compiler_3.6.2  magrittr_1.5    tools_3.6.2     htmltools_0.4.0
##  [5] yaml_2.2.0      Rcpp_1.0.2      stringi_1.4.3   rmarkdown_1.15 
##  [9] grid_3.6.2      stringr_1.4.0   xfun_0.10       digest_0.6.21  
## [13] rlang_0.4.0     lattice_0.20-38 evaluate_0.14