Wilke et al. richness of endemic diatoms
Three hundred and eighty diatom samples were analysed from the DEEP sediment succession, taken every 128 cm (successive samples are around 2,000–4,000 years apart) between 0 and 406.96 m c.d. and every 64 cm (around 2,000 years) between 406.96 and 446.65 m c.d. The higher sampling frequency in the early phase of Lake Ohrid was chosen because this period was environmentally very dynamic (Wagner et al. 2017, Panagiotopoulos et al. 2020). For comparison, we also analysed 123 samples from the phase before 1.36 Myr ago (446.65–583.92 m c.d.).
# Load required packages
library(coda)
library(foreach)
library(doParallel)
# Auxiliary functions
source("PyRate/pyrate_utilities.r")
source("AuxiliaryFunctions/lttFromMcmcLog.R")# Read data
Diatoms <- read.table("Data/Supplementary material_Diatoms.csv",
sep = "\t", header = TRUE, check.names = FALSE,
stringsAsFactors = FALSE)
# Exclude the endemic species we haven't found
Diatoms <- Diatoms[, -c(373:421)]
# Data cleaning
TaxonomyLifestyle <- Diatoms[1:12, -c(1:4)]
Diatoms <- Diatoms[-c(1:12), - c(1:3, 5)]
# Convert to numeric
Diatoms <- as.data.frame( apply(Diatoms, 2, as.numeric) )
# PyRate does not like special characters for species names:
SpeciesChangeNames <- data.frame(Species = colnames(Diatoms)[-1],
SpeciesPyRate = sprintf("SP_%03d", 1:(ncol(Diatoms)-1)),
stringsAsFactors = FALSE)
colnames(Diatoms) <- c("Age", SpeciesChangeNames$SpeciesPyRate)
TaxonomyLifestyle <- rbind(TaxonomyLifestyle, rep(NA, ncol(TaxonomyLifestyle)))
TaxonomyLifestyle[13, ] <- c("PyRateCode", SpeciesChangeNames$SpeciesPyRate)
# Endemic species
SpEndemic <- colnames(Diatoms)[-1][TaxonomyLifestyle[12, -1] == "Endemic"]
# Extinct endemic species
SpEx <- colnames(Diatoms)[-1][TaxonomyLifestyle[11, -1] == "Extinct" &
TaxonomyLifestyle[12, -1] == "Endemic"] Generate PyRate input
As described previously (Pires et al. 2018), the BDS analyses were informed on the endemic species that originated in the basin before the formation of Lake Ohrid 1.36 Myr ago and lack a robust age-control. They therefore did not factor in as speciation events, but the respective 11 extinction events were considered, as they took place 1.36–0 Myr ago.
SpeciationAfter136 <- rep(0, length(SpEndemic))
ExtinctionAfter136 <- SpeciationAfter136
Exclude <- SpeciationAfter136
for(i in 1:length(SpEndemic)){
Tmp <- Diatoms[, c("Age", SpEndemic[i])]
Observation <- which(Tmp[, 2] > 0)
FirstObservation <- Tmp$Age[max(Observation)]
LastObservation <- Tmp$Age[min(Observation)]
# Speciation between 1.36-0 Myr ago
if(!is.na(FirstObservation)){
SpeciationAfter136[i] <- 1
}
# Extinction between 1.36-0 Myr ago
if(!is.na(LastObservation) & SpEndemic[i] %in% SpEx){
ExtinctionAfter136[i] <- 1
}
# First and last observation before 1.36 Myr;
# no matter whether they are extant or extinct
# (Those species are not informative and do not contribute to the likelihood)
if(is.na(FirstObservation) & is.na(LastObservation)){
Exclude[i] <- 1
}
}
AliveBefore <- data.frame(taxon = SpEndemic,
speciation_in_bin = SpeciationAfter136,
extinction_in_bin = ExtinctionAfter136)
AliveBefore <- AliveBefore[Exclude == 0, ]
write.table(AliveBefore, "Results/BDS/AliveBefore.txt",
sep = "\t", row.names = FALSE, col.names = TRUE, quote = FALSE)We incorporated age uncertainties for the sediment succession (Wagner et al. 2019) in the inference of speciation and extinction times by generating 100 randomized inputs for the birth–death sampling (BDS) analyses with ages for the 380 diatom samples drawn from a uniform distribution of the mean age of the sample ± the age uncertainty of approximately 2,000 years.
Due to computational limits we only generate 2 randomized inputs for this reproducible example.
set.seed(23) # Reproducibility
DiaTrim <- Diatoms[, c("Age", SpEndemic[Exclude == 0])]
DiaTrim <- DiaTrim[-c(382:nrow(DiaTrim)), ]
# How many occurrences?
Rows <- sum(ifelse(DiaTrim[, -1] > 0, 1, 0) )
Counter <- 0
PyRateInput <- as.data.frame(matrix(NA, ncol = 5, nrow = Rows))
colnames(PyRateInput) <- c("Species", "Status", "MinT", "MaxT", "Trait")
PyRateInput$Trait <- 1
for(y in 2:ncol(DiaTrim)){ # Loop through all columns with species
for(i in 1:nrow(DiaTrim)){ # Loop through all rows
if(DiaTrim[i, y] > 0){ # If the species is present
Counter <- Counter + 1
PyRateInput[Counter, 1] <- colnames(DiaTrim)[y] # Species name
PyRateInput[Counter, 3] <- DiaTrim$Age[i] - 1.2
PyRateInput[Counter, 4] <- DiaTrim$Age[i] + 1.2
}
}
}
PyRateInput[, 2] <- "extant"
PyRateInput[PyRateInput$Species %in% SpEx, 2] <- "extinct"
# Rescale from 1400 kiloyears to 13 Million years
# for improving numerical stability and reasonable prior in PyRate
PyRateInput$MinT <- PyRateInput$MinT / 100
PyRateInput$MaxT <- PyRateInput$MaxT / 100
write.table(PyRateInput,
file = paste0("Results/BDS/TrimBounded.txt"),
sep = "\t", row.names = FALSE)
extract.ages(file = paste0("Results/BDS/TrimBounded.txt"), replicates = 2, outname = "")##
## replicate 1
## replicate 2
##
## PyRate input file was saved in: Results/BDS/TrimBounded.pyThe frequency of species’ occurrence in the current record may change over time due to ecosystem change altering the preservation potential or due to a more intensive screening during the early lake phase. Therefore, we defined three boundaries between four periods in which the mean sampling rate may differ according to major changes in δ18Olakewater and sampling intensity: 1.150, 0.866 and 0.455 Myr ago
PresShifts <- data.frame(Age = c(11.5, 8.66, 4.55))
write.table(PresShifts,
file = paste0("Results/BDS/PresShifts.txt"),
sep = "\t", col.names = FALSE, row.names = FALSE) Preservation models
Comparing preservation models (Silvestro et al. 2019) using the corrected Akaike information criterion (AICc) suggested a significant better fit to the fossil data of the time‑variable Poisson process for fossil sampling than a sampling rate that is either constant (ΔAICc = 72.1) or dependent on the lifespan of a species (ΔAICc = 425.2).
PresModelTest <- system(paste("python2 PyRate/PyRate.py",
"Results/BDS/TrimBounded.py",
"-PPmodeltest",
"-qShift Results/BDS/PresShifts.txt"),
intern = TRUE)
PresModelTest[grep("models", PresModelTest)]## [1] "models: ['HPP' 'NHPP' 'TPP']"PresModelTest[grep("AICc scores", PresModelTest)]## [1] "AICc scores: [ 3887.101 4240.858 3814.66 ]"Speciation and extinction times
We obtained posterior estimates of speciation and extinction times for each species using Markov Chain Monte Carlo (MCMC) sampling, which we ran for 5,000,000 iterations, sampling every 1,000 iterations and omitting the first 100,000 iterations as burn-in.
Due to computational limits we ran less iterations for this reproducible example.
# Define a call to python that we can run on several cores in parallel
runPyRate <- function(i){
system( paste("python2 PyRate/PyRate.py",
"Results/BDS/TrimBounded.py",
paste0("-j ", i),
"-n 500001", # Iterations to run
"-s 1000", # Sampling frequency
"-b 100000", # Burn-in
"-p 200000", # Print frequency
"-mG",
"-A 0",
"-mL 1",
"-mM 1",
"-trait_file Results/BDS/AliveBefore.txt", # see Pires et al. 2018
"-bound 13.66 0",
"-qShift Results/BDS/PresShifts.txt", # Preservation shifts
"-wd Results/BDS"))
}
# Run on 2 cores in parallel
Ncores <- 2
registerDoParallel(Ncores)
PyRate <- foreach(iter = 1:2, .inorder = FALSE) %dopar% runPyRate(iter)
stopImplicitCluster()Check effective sampling sizes
McmcLogFiles <- list.files("Results/BDS/pyrate_mcmc_logs", pattern = "GBD1-1_mcmc.log")
McmcLogs <- vector(mode = "list", length(McmcLogFiles))
for(i in 1:length(McmcLogFiles)){
McmcLogs[[i]] <- read.table(paste0("Results/BDS/pyrate_mcmc_logs/", McmcLogFiles[i]),
header = TRUE, sep = "\t")
}
lapply( McmcLogs, function(x) apply(x, 2, function(y) round(effectiveSize(as.mcmc(y)))))## [[1]]
## it posterior prior PP_lik BD_lik q_0
## 0 79 100 83 88 134
## q_1 q_2 q_3 alpha root_age death_age
## 187 111 181 261 31 0
## lambda_0 mu_0 tot_length SP_006_TS SP_008_TS SP_009_TS
## 401 724 27 31 113 173
## SP_010_TS SP_011_TS SP_017_TS SP_018_TS SP_020_TS SP_021_TS
## 401 31 31 269 139 31
## SP_022_TS SP_023_TS SP_024_TS SP_026_TS SP_028_TS SP_029_TS
## 401 82 31 31 31 9
## SP_031_TS SP_032_TS SP_033_TS SP_034_TS SP_035_TS SP_036_TS
## 14 140 107 31 10 341
## SP_038_TS SP_044_TS SP_045_TS SP_046_TS SP_047_TS SP_048_TS
## 165 265 326 31 379 192
## SP_051_TS SP_052_TS SP_053_TS SP_055_TS SP_056_TS SP_057_TS
## 418 31 251 66 66 122
## SP_060_TS SP_069_TS SP_070_TS SP_073_TS SP_074_TS SP_075_TS
## 31 8 335 31 83 31
## SP_076_TS SP_077_TS SP_078_TS SP_079_TS SP_081_TS SP_082_TS
## 31 31 31 259 31 273
## SP_085_TS SP_087_TS SP_089_TS SP_090_TS SP_092_TS SP_093_TS
## 16 31 401 401 31 31
## SP_094_TS SP_099_TS SP_101_TS SP_102_TS SP_104_TS SP_107_TS
## 266 162 277 184 31 87
## SP_108_TS SP_109_TS SP_115_TS SP_118_TS SP_121_TS SP_122_TS
## 57 186 68 108 146 315
## SP_125_TS SP_126_TS SP_127_TS SP_133_TS SP_134_TS SP_140_TS
## 342 118 101 31 31 273
## SP_145_TS SP_147_TS SP_148_TS SP_150_TS SP_156_TS SP_157_TS
## 188 570 226 222 33 54
## SP_162_TS SP_164_TS SP_166_TS SP_168_TS SP_170_TS SP_174_TS
## 156 31 31 31 246 401
## SP_177_TS SP_179_TS SP_181_TS SP_184_TS SP_185_TS SP_188_TS
## 31 286 31 216 31 31
## SP_189_TS SP_192_TS SP_193_TS SP_195_TS SP_196_TS SP_199_TS
## 171 31 30 269 111 32
## SP_207_TS SP_213_TS SP_214_TS SP_215_TS SP_216_TS SP_220_TS
## 31 401 93 174 294 31
## SP_222_TS SP_227_TS SP_231_TS SP_232_TS SP_233_TS SP_236_TS
## 401 31 31 175 317 401
## SP_242_TS SP_243_TS SP_247_TS SP_253_TS SP_262_TS SP_268_TS
## 401 332 155 146 31 31
## SP_269_TS SP_270_TS SP_272_TS SP_277_TS SP_280_TS SP_282_TS
## 31 304 31 274 8 401
## SP_285_TS SP_286_TS SP_287_TS SP_295_TS SP_302_TS SP_304_TS
## 401 31 357 31 31 21
## SP_305_TS SP_309_TS SP_310_TS SP_312_TS SP_313_TS SP_314_TS
## 31 201 125 333 31 31
## SP_318_TS SP_319_TS SP_320_TS SP_321_TS SP_332_TS SP_334_TS
## 31 131 31 31 31 31
## SP_338_TS SP_339_TS SP_340_TS SP_341_TS SP_349_TS SP_350_TS
## 31 49 29 259 328 220
## SP_351_TS SP_353_TS SP_354_TS SP_355_TS SP_357_TS SP_006_TE
## 212 401 314 179 235 260
## SP_008_TE SP_009_TE SP_010_TE SP_011_TE SP_017_TE SP_018_TE
## 140 0 0 0 0 0
## SP_020_TE SP_021_TE SP_022_TE SP_023_TE SP_024_TE SP_026_TE
## 0 0 0 0 0 0
## SP_028_TE SP_029_TE SP_031_TE SP_032_TE SP_033_TE SP_034_TE
## 0 0 0 0 0 0
## SP_035_TE SP_036_TE SP_038_TE SP_044_TE SP_045_TE SP_046_TE
## 0 333 401 0 0 0
## SP_047_TE SP_048_TE SP_051_TE SP_052_TE SP_053_TE SP_055_TE
## 0 0 0 0 0 0
## SP_056_TE SP_057_TE SP_060_TE SP_069_TE SP_070_TE SP_073_TE
## 0 0 0 0 18 401
## SP_074_TE SP_075_TE SP_076_TE SP_077_TE SP_078_TE SP_079_TE
## 0 401 137 91 0 401
## SP_081_TE SP_082_TE SP_085_TE SP_087_TE SP_089_TE SP_090_TE
## 401 0 0 0 0 0
## SP_092_TE SP_093_TE SP_094_TE SP_099_TE SP_101_TE SP_102_TE
## 0 0 0 0 280 140
## SP_104_TE SP_107_TE SP_108_TE SP_109_TE SP_115_TE SP_118_TE
## 0 401 0 0 0 0
## SP_121_TE SP_122_TE SP_125_TE SP_126_TE SP_127_TE SP_133_TE
## 0 0 201 0 0 0
## SP_134_TE SP_140_TE SP_145_TE SP_147_TE SP_148_TE SP_150_TE
## 0 0 0 0 0 28
## SP_156_TE SP_157_TE SP_162_TE SP_164_TE SP_166_TE SP_168_TE
## 0 0 0 0 0 0
## SP_170_TE SP_174_TE SP_177_TE SP_179_TE SP_181_TE SP_184_TE
## 401 0 0 0 0 0
## SP_185_TE SP_188_TE SP_189_TE SP_192_TE SP_193_TE SP_195_TE
## 0 0 0 0 122 0
## SP_196_TE SP_199_TE SP_207_TE SP_213_TE SP_214_TE SP_215_TE
## 0 0 401 0 0 0
## SP_216_TE SP_220_TE SP_222_TE SP_227_TE SP_231_TE SP_232_TE
## 0 0 0 0 0 310
## SP_233_TE SP_236_TE SP_242_TE SP_243_TE SP_247_TE SP_253_TE
## 145 0 0 0 0 0
## SP_262_TE SP_268_TE SP_269_TE SP_270_TE SP_272_TE SP_277_TE
## 191 401 401 0 0 0
## SP_280_TE SP_282_TE SP_285_TE SP_286_TE SP_287_TE SP_295_TE
## 0 0 0 0 0 0
## SP_302_TE SP_304_TE SP_305_TE SP_309_TE SP_310_TE SP_312_TE
## 0 0 0 0 0 0
## SP_313_TE SP_314_TE SP_318_TE SP_319_TE SP_320_TE SP_321_TE
## 0 0 0 0 0 0
## SP_332_TE SP_334_TE SP_338_TE SP_339_TE SP_340_TE SP_341_TE
## 0 0 401 157 61 204
## SP_349_TE SP_350_TE SP_351_TE SP_353_TE SP_354_TE SP_355_TE
## 0 0 0 39 0 0
## SP_357_TE
## 18
##
## [[2]]
## it posterior prior PP_lik BD_lik q_0
## 0 90 97 99 34 98
## q_1 q_2 q_3 alpha root_age death_age
## 161 145 165 186 56 0
## lambda_0 mu_0 tot_length SP_006_TS SP_008_TS SP_009_TS
## 401 401 22 56 249 87
## SP_010_TS SP_011_TS SP_017_TS SP_018_TS SP_020_TS SP_021_TS
## 401 56 56 263 185 56
## SP_022_TS SP_023_TS SP_024_TS SP_026_TS SP_028_TS SP_029_TS
## 335 57 56 56 56 8
## SP_031_TS SP_032_TS SP_033_TS SP_034_TS SP_035_TS SP_036_TS
## 15 73 81 56 7 315
## SP_038_TS SP_044_TS SP_045_TS SP_046_TS SP_047_TS SP_048_TS
## 235 224 344 56 340 246
## SP_051_TS SP_052_TS SP_053_TS SP_055_TS SP_056_TS SP_057_TS
## 401 56 239 85 43 158
## SP_060_TS SP_069_TS SP_070_TS SP_073_TS SP_074_TS SP_075_TS
## 56 10 288 56 128 56
## SP_076_TS SP_077_TS SP_078_TS SP_079_TS SP_081_TS SP_082_TS
## 56 56 56 326 56 287
## SP_085_TS SP_087_TS SP_089_TS SP_090_TS SP_092_TS SP_093_TS
## 12 56 323 401 56 56
## SP_094_TS SP_099_TS SP_101_TS SP_102_TS SP_104_TS SP_107_TS
## 334 101 231 126 56 132
## SP_108_TS SP_109_TS SP_115_TS SP_118_TS SP_121_TS SP_122_TS
## 161 151 117 110 106 297
## SP_125_TS SP_126_TS SP_127_TS SP_133_TS SP_134_TS SP_140_TS
## 401 142 124 56 56 385
## SP_145_TS SP_147_TS SP_148_TS SP_150_TS SP_156_TS SP_157_TS
## 158 736 254 181 97 93
## SP_162_TS SP_164_TS SP_166_TS SP_168_TS SP_170_TS SP_174_TS
## 82 56 56 56 254 401
## SP_177_TS SP_179_TS SP_181_TS SP_184_TS SP_185_TS SP_188_TS
## 56 284 56 232 56 56
## SP_189_TS SP_192_TS SP_193_TS SP_195_TS SP_196_TS SP_199_TS
## 228 56 34 234 160 3
## SP_207_TS SP_213_TS SP_214_TS SP_215_TS SP_216_TS SP_220_TS
## 56 309 48 268 245 56
## SP_222_TS SP_227_TS SP_231_TS SP_232_TS SP_233_TS SP_236_TS
## 401 56 56 162 401 401
## SP_242_TS SP_243_TS SP_247_TS SP_253_TS SP_262_TS SP_268_TS
## 401 401 167 133 56 56
## SP_269_TS SP_270_TS SP_272_TS SP_277_TS SP_280_TS SP_282_TS
## 56 215 56 298 13 401
## SP_285_TS SP_286_TS SP_287_TS SP_295_TS SP_302_TS SP_304_TS
## 401 56 218 56 56 5
## SP_305_TS SP_309_TS SP_310_TS SP_312_TS SP_313_TS SP_314_TS
## 56 253 114 338 56 56
## SP_318_TS SP_319_TS SP_320_TS SP_321_TS SP_332_TS SP_334_TS
## 56 200 56 56 56 56
## SP_338_TS SP_339_TS SP_340_TS SP_341_TS SP_349_TS SP_350_TS
## 56 97 103 189 277 299
## SP_351_TS SP_353_TS SP_354_TS SP_355_TS SP_357_TS SP_006_TE
## 172 401 401 141 265 338
## SP_008_TE SP_009_TE SP_010_TE SP_011_TE SP_017_TE SP_018_TE
## 283 0 0 0 0 0
## SP_020_TE SP_021_TE SP_022_TE SP_023_TE SP_024_TE SP_026_TE
## 0 0 0 0 0 0
## SP_028_TE SP_029_TE SP_031_TE SP_032_TE SP_033_TE SP_034_TE
## 0 0 0 0 0 0
## SP_035_TE SP_036_TE SP_038_TE SP_044_TE SP_045_TE SP_046_TE
## 0 401 295 0 0 0
## SP_047_TE SP_048_TE SP_051_TE SP_052_TE SP_053_TE SP_055_TE
## 0 0 0 0 0 0
## SP_056_TE SP_057_TE SP_060_TE SP_069_TE SP_070_TE SP_073_TE
## 0 0 0 0 32 401
## SP_074_TE SP_075_TE SP_076_TE SP_077_TE SP_078_TE SP_079_TE
## 0 337 140 132 0 401
## SP_081_TE SP_082_TE SP_085_TE SP_087_TE SP_089_TE SP_090_TE
## 401 0 0 0 0 0
## SP_092_TE SP_093_TE SP_094_TE SP_099_TE SP_101_TE SP_102_TE
## 0 0 0 0 149 401
## SP_104_TE SP_107_TE SP_108_TE SP_109_TE SP_115_TE SP_118_TE
## 0 334 0 0 0 0
## SP_121_TE SP_122_TE SP_125_TE SP_126_TE SP_127_TE SP_133_TE
## 0 0 129 0 0 0
## SP_134_TE SP_140_TE SP_145_TE SP_147_TE SP_148_TE SP_150_TE
## 0 0 0 0 0 36
## SP_156_TE SP_157_TE SP_162_TE SP_164_TE SP_166_TE SP_168_TE
## 0 0 0 0 0 0
## SP_170_TE SP_174_TE SP_177_TE SP_179_TE SP_181_TE SP_184_TE
## 401 0 0 0 0 0
## SP_185_TE SP_188_TE SP_189_TE SP_192_TE SP_193_TE SP_195_TE
## 0 0 0 0 158 0
## SP_196_TE SP_199_TE SP_207_TE SP_213_TE SP_214_TE SP_215_TE
## 0 0 326 0 0 0
## SP_216_TE SP_220_TE SP_222_TE SP_227_TE SP_231_TE SP_232_TE
## 0 0 0 0 0 291
## SP_233_TE SP_236_TE SP_242_TE SP_243_TE SP_247_TE SP_253_TE
## 135 0 0 0 0 0
## SP_262_TE SP_268_TE SP_269_TE SP_270_TE SP_272_TE SP_277_TE
## 456 401 401 0 0 0
## SP_280_TE SP_282_TE SP_285_TE SP_286_TE SP_287_TE SP_295_TE
## 0 0 0 0 0 0
## SP_302_TE SP_304_TE SP_305_TE SP_309_TE SP_310_TE SP_312_TE
## 0 0 0 0 0 0
## SP_313_TE SP_314_TE SP_318_TE SP_319_TE SP_320_TE SP_321_TE
## 0 0 0 0 0 0
## SP_332_TE SP_334_TE SP_338_TE SP_339_TE SP_340_TE SP_341_TE
## 0 0 470 152 51 137
## SP_349_TE SP_350_TE SP_351_TE SP_353_TE SP_354_TE SP_355_TE
## 0 0 0 82 0 0
## SP_357_TE
## 10Endemic richness through time
The accumulation of endemic species richness through time is the cumulative sum of all origination minus the cumulative sum of all extinction times in a single MCMC sample.
ColTs <- range(which(grepl("TS", colnames(McmcLogs[[1]]))))
ColTs <- ColTs[1]:ColTs[2]
ColTe <- range(which(grepl("TE", colnames(McmcLogs[[1]]))))
ColTe <- ColTe[1]:ColTe[2]
Ltt <- vector(mode = "list", length(McmcLogFiles))
for(i in 1:length(McmcLogs)){
Ltt[[i]] <- apply(McmcLogs[[i]], 1, function(x) lttFromMcmcLog(x, ColTs, ColTe))
}
Ltt <- do.call(cbind, Ltt)
Keep <- apply(Ltt, 2, function(x) !all(is.na(x))) # Omit columns with all NAs
Ltt <- Ltt[, Keep]Endemic richness was averaged over all 100 BDS analyses and the uncertainty summarized by highest posterior density intervals with probabilities masses of 0.5−0.95 incremented by 0.05.
Prob <- seq(0.95, 0.05, by = -0.05)
LttHpd <- vector(mode = "list", length(Prob))
for(i in 1:length(Prob)){
LttHpd[[i]] <- apply(Ltt[-nrow(Ltt), ], 1, function(x) HPDinterval(coda:::mcmc(x), prob = Prob[i]))
}
MeanDiv <- rowMeans(Ltt[-nrow(Ltt), ])
TimeSeq <- seq(0, 13.7, 0.05) / 10
par(las = 1, mar = c(5.1, 4.5, 0.1, 0.8))
plot(seq(0, 14, 2)/10, rep(1, 8),
type = "n", xaxs = "i",
xlim = c(1.4, 0), ylim = c(50, 130),
xlab = "Age (Myr ago)", ylab = "Endemic species richness (N)")
for(i in 1:length(Prob)){
polygon(c(TimeSeq[-length(TimeSeq)], rev(TimeSeq[-length(TimeSeq)])),
c(LttHpd[[i]][1, ], rev(LttHpd[[i]][2, ])),
col = adjustcolor("green", alpha = 1/length(Prob)), border = NA)
}
lines(TimeSeq[-length(TimeSeq)], MeanDiv, col = "green")
Plattform and resources
Sys.time()## [1] "2019-12-31 20:31:05 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] parallel stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] doParallel_1.0.15 iterators_1.0.12 foreach_1.4.7 coda_0.19-3
##
## loaded via a namespace (and not attached):
## [1] Rcpp_1.0.2 codetools_0.2-16 lattice_0.20-38 digest_0.6.21
## [5] grid_3.6.2 magrittr_1.5 evaluate_0.14 rlang_0.4.0
## [9] stringi_1.4.3 rmarkdown_1.15 tools_3.6.2 stringr_1.4.0
## [13] xfun_0.10 yaml_2.2.0 compiler_3.6.2 htmltools_0.4.0
## [17] knitr_1.25