hitAPI <- FALSE modifiedDate <- '2018.07.31' # download species measurement from PBDB (or load previous download) if(hitAPI) { rawData <- read.csv(file=paste("https://paleobiodb.org/data1.2/specs/measurements.csv?base_name=Brachiopoda&taxon_reso=species&interval=Permian&show=spec,class,ref,entname&specs_modified_before=", modifiedDate, sep=""), header=TRUE) allGenera <- paste(levels(rawData$genus), collapse=",") genusData <- read.csv(file=paste("https://paleobiodb.org/data1.2/occs/list.csv?base_name=Brachiopoda&idreso=species&show=class&occs_modified_before=", modifiedDate, sep=""), header=TRUE) genusData <- subset(genusData, is.element(genus, levels(rawData$genus))) genusData$accepted_name <- factor(as.character(genusData$accepted_name)) genusData$genus <- factor(as.character(genusData$genus)) genusRanges <- data.frame('genus'=levels(genusData$genus), 'firstapp_max_ma'=tapply(genusData$max_ma, genusData$genus, max, na.rm=TRUE), 'lastapp_min_ma'=tapply(genusData$min_ma, genusData$genus, min, na.rm=TRUE), 'n_species'=tapply(genusData$accepted_name, genusData$genus, function(x){return(length(unique(as.character(x))))}) ) } else { load('permianBrachs.RData') } # created data frame of species measurement and specimen ages speciesSizes <- data.frame('genus'=factor(tapply(rawData$genus, rawData$accepted_name, function(x){return(as.character(x[1]))})), 'maxAge'=tapply(rawData$max_ma, rawData$accepted_name, max, na.rm=TRUE),'minAge'=tapply(rawData$min_ma, rawData$accepted_name, min, na.rm=TRUE), tapply(rawData$average, list(rawData$accepted_name, rawData$measurement), max, na.rm=TRUE), 'maxDim' = tapply(rawData$average, rawData$accepted_name, max, na.rm=TRUE)) # create data frame genusSizes <- data.frame('nSpecies'=as.numeric(table(speciesSizes$genus)), 'maxAge'=tapply(speciesSizes$maxAge, speciesSizes$genus, max, na.rm=TRUE), 'minAge'=tapply(speciesSizes$minAge, speciesSizes$genus, min, na.rm=TRUE), 'duration'=tapply(speciesSizes$maxAge, speciesSizes$genus, max, na.rm=TRUE) - tapply(speciesSizes$minAge, speciesSizes$genus, min, na.rm=TRUE), 'meanSize'=tapply(speciesSizes$maxDim, speciesSizes$genus, mean), 'minSize'=tapply(speciesSizes$maxDim, speciesSizes$genus, min), 'maxSize'=tapply(speciesSizes$maxDim, speciesSizes$genus, max) ) # ranges from speciesSizes are just for the Permian occurrences, these are for all genus occurrences genusSizes$genusFad <- genusRanges$firstapp_max_ma[match(genusRanges$genus, rownames(genusSizes))] genusSizes$genusLad <- genusRanges$lastapp_min_ma[match(genusRanges$genus, rownames(genusSizes))] genusSizes$genusDuration <- genusSizes$genusFad - genusSizes$genusLad genusSizes$totalSpecies <- genusRanges$n_species[match(genusRanges$genus, rownames(genusSizes))] setEPS() postscript(file='Fig6_Payne_2col.eps', height=5.8, width=5.8) par(mfrow=c(2,2), mar=c(4,3.75,2,1), pch=16, las=1, mgp=c(2.5,1,0)) cexAxis <- 1 cexLab <- 1 cexMain <- 1 cexText <- 1 cexStat <- 1 yMinBD <- 0.2 mySep <- ", " # based only on Permian species with sizes plot(log2(genusSizes$nSpecies), log10(genusSizes$meanSize), xlab=expression(paste("log"[2]," Number of species",sep="")), ylab=expression(paste("Body size (log"[10]," mm)", sep="")), main="Species per Genus (Permian)", cex.axis=cexAxis, cex.lab=cexLab, cex.main=cexMain) mtext("A", side=3, font=2, cex=cexText, line=0.75, adj=0, at=-2) myCor <- cor.test(log2(genusSizes$nSpecies), log10(genusSizes$meanSize), method="spearman") text(max(log2(genusSizes$nSpecies)), min(log10(genusSizes$meanSize)), labels=bquote(rho == .(round(myCor$estimate,4)) ~ .(mySep) ~ 'p' == .(round(myCor$p.value,4))), cex=cexStat <- 1, adj=1) plot(log10(genusSizes$duration), log10(genusSizes$meanSize), ylim=c(yMinBD, max(log10(genusSizes$meanSize))), xlab=expression(paste("Genus duration (log"[10]," Ma)", sep="")), ylab=expression(paste("Body size (log"[10]," mm)", sep="")), main="Duration (Permian)", cex.axis=cexAxis, cex.lab=cexLab, cex.main=cexMain) mtext("B", side=3, font=2, cex=cexText, line=0.75, adj=0, at=-0.15) myCor <- cor.test(log10(genusSizes$duration), log10(genusSizes$meanSize), method="spearman") text(max(log10(genusSizes$duration)), yMinBD, labels=bquote(rho == .(round(myCor$estimate,4)) ~ .(mySep) ~ 'p' == .(round(myCor$p.value,4))), cex=cexStat <- 1, adj=1) # based on total number of species and duration plot(log2(genusSizes$totalSpecies), log10(genusSizes$meanSize), xlab=expression(paste("log"[2]," Number of species",sep="")), ylab=expression(paste("Body size (log"[10]," mm)", sep="")), main="Species per Genus (total)", cex.axis=cexAxis, cex.lab=cexLab, cex.main=cexMain) mtext("C", side=3, font=2, cex=cexText, line=0.75, adj=0, at=-2) myCor <- cor.test(log2(genusSizes$totalSpecies), log10(genusSizes$meanSize), method="spearman") text(max(log2(genusSizes$totalSpecies)), min(log10(genusSizes$meanSize)), labels=bquote(rho == .(round(myCor$estimate,4)) ~ .(mySep) ~ 'p' == .(round(myCor$p.value,4))), cex=cexStat <- 1, adj=1) plot(log10(genusSizes$genusDuration), log10(genusSizes$meanSize), ylim=c(yMinBD, max(log10(genusSizes$meanSize))), xlab=expression(paste("Genus duration (log"[10]," Ma)", sep="")), ylab=expression(paste("Body size (log"[10]," mm)", sep="")), main="Duration (total)", cex.axis=cexAxis, cex.lab=cexLab, cex.main=cexMain) mtext("D", side=3, font=2, cex=cexText, line=0.75, adj=0, at=-0.47) myCor <- cor.test(log10(genusSizes$genusDuration), log10(genusSizes$meanSize), method="spearman") text(max(log10(genusSizes$genusDuration)), yMinBD, labels=bquote(rho == .(round(myCor$estimate,4)) ~ .(mySep) ~ 'p' == .(round(myCor$p.value,4))), cex=cexStat <- 1, adj=1) dev.off() summary(glm(log10(meanSize) ~ log2(nSpecies) + log10(duration), data=genusSizes)) summary(glm(log10(meanSize) ~ log2(totalSpecies) + log10(genusDuration), data=genusSizes)) save(rawData, speciesSizes, genusSizes, genusRanges, file='permianBrachs.RData')