####10,000 iterations of a null model by conducting a random permutation of ecological trends by shuffling regions then followed by Spearman’s rank correlation tests #(note that this script takes at least 26 hours to run) ##Loop, plots, and Spearman's tests of ecological percentages of PBDB Jurassic regional occurrences #Read percent confidence interval function source("percentConfidenceInterval.R") #Read percent ecology function to calculate percentages with confidence intervals source("percentEcolTrait.R") #Read dataset occurrences2 <- read.csv('ClusteredOccurrences.csv', header=TRUE, stringsAsFactors=FALSE) #Create object to hold results of 10,000 iteration null model hold <- vector(mode = "numeric", length = 10000) EU_hold <- vector(mode = "numeric", length = 10000) #Create objects to define dimensions of matrices below stagecolumn <- 4 JurassicStages <- sort(unique(occurrences2[occurrences2$stage!="0", stagecolumn])) geoClusters <- sort(unique(occurrences2$geo_cluster)) #Loop to calculate ecological trends by stage per region for (i in 1:10000) { geoClusters <- sample(geoClusters, length(geoClusters)) #regional shuffle #Create empty matrices that will be used to fill estimated percentages calculated by running loop sinfaunalTableE <- matrix(data=0, nrow=length(JurassicStages), ncol=length(geoClusters), dimnames=list(JurassicStages, geoClusters)) dinfaunalTableE <- matrix(data=0, nrow=length(JurassicStages), ncol=length(geoClusters), dimnames=list(JurassicStages, geoClusters)) epifaunalTableE <- matrix(data=0, nrow=length(JurassicStages), ncol=length(geoClusters), dimnames=list(JurassicStages, geoClusters)) cementedTableE <- matrix(data=0, nrow=length(JurassicStages), ncol=length(geoClusters), dimnames=list(JurassicStages, geoClusters)) freelivingTableE <- matrix(data=0, nrow=length(JurassicStages), ncol=length(geoClusters), dimnames=list(JurassicStages, geoClusters)) pedunculateTableE <- matrix(data=0, nrow=length(JurassicStages), ncol=length(geoClusters), dimnames=list(JurassicStages, geoClusters)) mobileTableE <- matrix(data=0, nrow=length(JurassicStages), ncol=length(geoClusters), dimnames=list(JurassicStages, geoClusters)) limitedmobileTableE <- matrix(data=0, nrow=length(JurassicStages), ncol=length(geoClusters), dimnames=list(JurassicStages, geoClusters)) sessileTableE <- matrix(data=0, nrow=length(JurassicStages), ncol=length(geoClusters), dimnames=list(JurassicStages, geoClusters)) #Create empty matrices that will be used to fill lower confidence limit percentages calculated by running loop sinfaunalTableLCL <- matrix(data=0, nrow=length(JurassicStages), ncol=length(geoClusters), dimnames=list(JurassicStages, geoClusters)) dinfaunalTableLCL <- matrix(data=0, nrow=length(JurassicStages), ncol=length(geoClusters), dimnames=list(JurassicStages, geoClusters)) epifaunalTableLCL <- matrix(data=0, nrow=length(JurassicStages), ncol=length(geoClusters), dimnames=list(JurassicStages, geoClusters)) cementedTableLCL <- matrix(data=0, nrow=length(JurassicStages), ncol=length(geoClusters), dimnames=list(JurassicStages, geoClusters)) freelivingTableLCL <- matrix(data=0, nrow=length(JurassicStages), ncol=length(geoClusters), dimnames=list(JurassicStages, geoClusters)) pedunculateTableLCL <- matrix(data=0, nrow=length(JurassicStages), ncol=length(geoClusters), dimnames=list(JurassicStages, geoClusters)) mobileTableLCL <- matrix(data=0, nrow=length(JurassicStages), ncol=length(geoClusters), dimnames=list(JurassicStages, geoClusters)) limitedmobileTableLCL <- matrix(data=0, nrow=length(JurassicStages), ncol=length(geoClusters), dimnames=list(JurassicStages, geoClusters)) sessileTableLCL <- matrix(data=0, nrow=length(JurassicStages), ncol=length(geoClusters), dimnames=list(JurassicStages, geoClusters)) #Create empty matrices that will be used to fill upper condfidence limit percentages calculated by running loop sinfaunalTableUCL <- matrix(data=0, nrow=length(JurassicStages), ncol=length(geoClusters), dimnames=list(JurassicStages, geoClusters)) dinfaunalTableUCL <- matrix(data=0, nrow=length(JurassicStages), ncol=length(geoClusters), dimnames=list(JurassicStages, geoClusters)) epifaunalTableUCL <- matrix(data=0, nrow=length(JurassicStages), ncol=length(geoClusters), dimnames=list(JurassicStages, geoClusters)) cementedTableUCL <- matrix(data=0, nrow=length(JurassicStages), ncol=length(geoClusters), dimnames=list(JurassicStages, geoClusters)) freelivingTableUCL <- matrix(data=0, nrow=length(JurassicStages), ncol=length(geoClusters), dimnames=list(JurassicStages, geoClusters)) pedunculateTableUCL <- matrix(data=0, nrow=length(JurassicStages), ncol=length(geoClusters), dimnames=list(JurassicStages, geoClusters)) mobileTableUCL <- matrix(data=0, nrow=length(JurassicStages), ncol=length(geoClusters), dimnames=list(JurassicStages, geoClusters)) limitedmobileTableUCL <- matrix(data=0, nrow=length(JurassicStages), ncol=length(geoClusters), dimnames=list(JurassicStages, geoClusters)) sessileTableUCL <- matrix(data=0, nrow=length(JurassicStages), ncol=length(geoClusters), dimnames=list(JurassicStages, geoClusters)) #Loop to calculate percentage of each ecological trait by cluster for each stage for (time in 1:length(JurassicStages)) { for (geo in 1:length(geoClusters)) { sinfaunal <- percentEcolTrait(occurrences2, "life_habit", "shallow infaunal", desiredStages=JurassicStages[time], desiredGeoclusters=geoClusters[geo]) sinfaunalTableE[time, geo] <- sinfaunal[2] sinfaunalTableLCL[time, geo] <- sinfaunal[1] sinfaunalTableUCL[time, geo] <- sinfaunal[3] dinfaunal <- percentEcolTrait(occurrences2, "life_habit", "deep infaunal", desiredStages=JurassicStages[time], desiredGeoclusters=geoClusters[geo]) dinfaunalTableE[time, geo] <- dinfaunal[2] dinfaunalTableLCL[time, geo] <- dinfaunal[1] dinfaunalTableUCL[time, geo] <- dinfaunal[3] epifaunal <- percentEcolTrait(occurrences2, "life_habit", "epifaunal", desiredStages=JurassicStages[time], desiredGeoclusters=geoClusters[geo]) epifaunalTableE[time, geo] <- epifaunal[2] epifaunalTableLCL[time, geo] <- epifaunal[1] epifaunalTableUCL[time, geo] <- epifaunal[3] cemented <- percentEcolTrait(occurrences2, "stability", "cemented", desiredStages=JurassicStages[time], desiredGeoclusters=geoClusters[geo]) cementedTableE[time, geo] <- cemented[2] cementedTableLCL[time, geo] <- cemented[1] cementedTableUCL[time, geo] <- cemented[3] freeliving <- percentEcolTrait(occurrences2, "stability", "free-living", desiredStages=JurassicStages[time], desiredGeoclusters=geoClusters[geo]) freelivingTableE[time, geo] <- freeliving[2] freelivingTableLCL[time, geo] <- freeliving[1] freelivingTableUCL[time, geo] <- freeliving[3] pedunculate <- percentEcolTrait(occurrences2, "stability", "pedunculate", desiredStages=JurassicStages[time], desiredGeoclusters=geoClusters[geo]) pedunculateTableE[time, geo] <- pedunculate[2] pedunculateTableLCL[time, geo] <- pedunculate[1] pedunculateTableUCL[time, geo] <- pedunculate[3] mobile <- percentEcolTrait(occurrences2, "motility", "mobile", desiredStages=JurassicStages[time], desiredGeoclusters=geoClusters[geo]) mobileTableE[time, geo] <- mobile[2] mobileTableLCL[time, geo] <- mobile[1] mobileTableUCL[time, geo] <- mobile[3] limitedmobile <- percentEcolTrait(occurrences2, "motility", "limited mobile", desiredStages=JurassicStages[time], desiredGeoclusters=geoClusters[geo]) limitedmobileTableE[time, geo] <- limitedmobile[2] limitedmobileTableLCL[time, geo] <- limitedmobile[1] limitedmobileTableUCL[time, geo] <- limitedmobile[3] sessile <- percentEcolTrait(occurrences2, "motility", "sessile", desiredStages=JurassicStages[time], desiredGeoclusters=geoClusters[geo]) sessileTableE[time, geo] <- sessile[2] sessileTableLCL[time, geo] <- sessile[1] sessileTableUCL[time, geo] <- sessile[3] } } #Create vector of Jurassic stage boundary ages (Myr) based on the 2015 ICS Chronostratigraphy Chart stageBoundaryAges <- c(201.3, 199.3, 190.8, 182.7, 174.1, 170.3, 168.3, 166.1, 163.5, 157.3, 152.1, 145.0) #Duration of each stage (Myr) – Hettangian: 2.0, Sinemurian: 8.5, Pliensbachian: 8.1, Toarcian: 8.6, Aalenian: 3.8, Bajocian: 2.0, Bathonian: 2.2, Callovian: 2.6, Oxfordian: 6.2, Kimmeridian: 5.2, Tithonian: 7.1 #Calculate midpoint of each stage to include in plots midpointAges <- stageBoundaryAges[-length(stageBoundaryAges)] + diff(stageBoundaryAges)/2 midpointAges <- sort(midpointAges) #Create vector to hold number each p-value pvalues <- vector(mode="numeric", length=45) EU_pvalues <- vector(mode="numeric", length=9) #Spearman's rank correlation tests for Europe Cluster test1 <- cor.test(midpointAges, sinfaunalTableE[,1], method='spearman') test2 <- cor.test(midpointAges, dinfaunalTableE[,1], method='spearman') test3 <- cor.test(midpointAges, epifaunalTableE[,1], method='spearman') test4 <- cor.test(midpointAges, freelivingTableE[,1], method='spearman') test5 <- cor.test(midpointAges, cementedTableE[,1], method='spearman') test6 <- cor.test(midpointAges, pedunculateTableE[,1], method='spearman') test7 <- cor.test(midpointAges, sessileTableE[,1], method='spearman') test8 <- cor.test(midpointAges, limitedmobileTableE[,1], method='spearman') test9 <- cor.test(midpointAges, mobileTableE[,1], method='spearman') #Spearman's rank correlation tests for Middle East Cluster test10 <- cor.test(midpointAges, sinfaunalTableE[,2], method='spearman') test11 <- cor.test(midpointAges, dinfaunalTableE[,2], method='spearman') test12 <- cor.test(midpointAges, epifaunalTableE[,2], method='spearman') test13 <- cor.test(midpointAges, freelivingTableE[,2], method='spearman') test14 <- cor.test(midpointAges, cementedTableE[,2], method='spearman') test15 <- cor.test(midpointAges, pedunculateTableE[,2], method='spearman') test16 <- cor.test(midpointAges, sessileTableE[,2], method='spearman') test17 <- cor.test(midpointAges, limitedmobileTableE[,2], method='spearman') test18 <- cor.test(midpointAges, mobileTableE[,2], method='spearman') #Spearman's rank correlation tests for New Zealand Cluster test19 <- cor.test(midpointAges, sinfaunalTableE[,3], method='spearman') test20 <- cor.test(midpointAges, dinfaunalTableE[,3], method='spearman') test21 <- cor.test(midpointAges, epifaunalTableE[,3], method='spearman') test22 <- cor.test(midpointAges, freelivingTableE[,3], method='spearman') test23 <- cor.test(midpointAges, cementedTableE[,3], method='spearman') test24 <- cor.test(midpointAges, pedunculateTableE[,3], method='spearman') test25 <- cor.test(midpointAges, sessileTableE[,3], method='spearman') test26 <- cor.test(midpointAges, limitedmobileTableE[,3], method='spearman') test27 <- cor.test(midpointAges, mobileTableE[,3], method='spearman') #Spearman's rank correlation tests for North America Cluster test28 <- cor.test(midpointAges, sinfaunalTableE[,4], method='spearman') test29 <- cor.test(midpointAges, dinfaunalTableE[,4], method='spearman') test30 <- cor.test(midpointAges, epifaunalTableE[,4], method='spearman') test31 <- cor.test(midpointAges, freelivingTableE[,4], method='spearman') test32 <- cor.test(midpointAges, cementedTableE[,4], method='spearman') test33 <- cor.test(midpointAges, pedunculateTableE[,4], method='spearman') test34 <- cor.test(midpointAges, sessileTableE[,4], method='spearman') test35 <- cor.test(midpointAges, limitedmobileTableE[,4], method='spearman') test36 <- cor.test(midpointAges, mobileTableE[,4], method='spearman') #Spearman's rank correlation tests for South America Cluster test37 <- cor.test(midpointAges, sinfaunalTableE[,5], method='spearman') test38 <- cor.test(midpointAges, dinfaunalTableE[,5], method='spearman') test39 <- cor.test(midpointAges, epifaunalTableE[,5], method='spearman') test40 <- cor.test(midpointAges, freelivingTableE[,5], method='spearman') test41 <- cor.test(midpointAges, cementedTableE[,5], method='spearman') test42 <- cor.test(midpointAges, pedunculateTableE[,5], method='spearman') test43 <- cor.test(midpointAges, sessileTableE[,5], method='spearman') test44 <- cor.test(midpointAges, limitedmobileTableE[,5], method='spearman') test45 <- cor.test(midpointAges, mobileTableE[,5], method='spearman') #Fill pvalues vector with p-values of each Speaman's correlation test pvalues[1] <- test1$p.value pvalues[2] <- test2$p.value pvalues[3] <- test3$p.value pvalues[4] <- test4$p.value pvalues[5] <- test5$p.value pvalues[6] <- test6$p.value pvalues[7] <- test7$p.value pvalues[8] <- test8$p.value pvalues[9] <- test9$p.value pvalues[10] <- test10$p.value pvalues[11] <- test11$p.value pvalues[12] <- test12$p.value pvalues[13] <- test13$p.value pvalues[14] <- test14$p.value pvalues[15] <- test15$p.value pvalues[16] <- test16$p.value pvalues[17] <- test17$p.value pvalues[18] <- test18$p.value pvalues[19] <- test19$p.value pvalues[20] <- test20$p.value pvalues[21] <- test21$p.value pvalues[22] <- test22$p.value pvalues[23] <- test23$p.value pvalues[24] <- test24$p.value pvalues[25] <- test25$p.value pvalues[26] <- test26$p.value pvalues[27] <- test27$p.value pvalues[28] <- test28$p.value pvalues[29] <- test29$p.value pvalues[30] <- test30$p.value pvalues[31] <- test31$p.value pvalues[32] <- test32$p.value pvalues[33] <- test33$p.value pvalues[34] <- test34$p.value pvalues[35] <- test35$p.value pvalues[36] <- test36$p.value pvalues[37] <- test37$p.value pvalues[38] <- test38$p.value pvalues[39] <- test39$p.value pvalues[40] <- test40$p.value pvalues[41] <- test41$p.value pvalues[42] <- test42$p.value pvalues[43] <- test43$p.value pvalues[44] <- test44$p.value pvalues[45] <- test45$p.value #Fill vector with number of pvalues at or lesser than 0.05 for each iteration hold[i] <- length(pvalues[pvalues <= 0.05]) EU_pvalues[1] <- test1$p.value EU_pvalues[2] <- test2$p.value EU_pvalues[3] <- test3$p.value EU_pvalues[4] <- test4$p.value EU_pvalues[5] <- test5$p.value EU_pvalues[6] <- test6$p.value EU_pvalues[7] <- test7$p.value EU_pvalues[8] <- test8$p.value EU_pvalues[9] <- test9$p.value EU_hold[i] <- length(EU_pvalues[EU_pvalues <= 0.05]) } #Create file with histogram of randomized iterations pdf("histogramGlobal.pdf", width=7, height=7) hist(hold, xlab="# of sig. trends", main="10000 Trials (Global)", las=1, col="gray") abline(v=8, col="red") dev.off() pdf("histogramEurope.pdf", width=7, height=7) hist(EU_hold, xlab="# of sig. trends", main="10000 Trials (Europe)", las=1, col="gray") abline(v=4, col="red") dev.off() #Calculate p-values length(hold[hold >= 8])/10000 length(EU_hold[EU_hold >= 4])/10000