
# Code for plotting the figures in Fritz et al.

# Main Figures
# 1 (line   30): sampling maps for mammalian locations (in grid cells) and
#                paleobotanical locations (only those used to calculate NPP)
# 2 (line  230): temporal diversity patterns for continents and regions
#                and temporal NPP patterns for continents (line 435)
# 3 (line  520): diversity-climate models and relationships
#                (scatterplots of fitted values) and comparison to present day

# Supplementary Figures
# S1 (line  670): raw genus counts, range-through genera, and the number of locations
# S2 (line  890): correlation of diversity estimators, and scatterplots of
#                 observed + range-through genera vs gamma diversity estimates
# S3 (line 1035): gamma diversity patterns with different location definitions
# S4 (line 1200): gamma estimation with the present day data (assessing
#                 performance where true richness is known)


# S5 was done in Illustrator (conceptual illustration)
# S6 and S7 were done by Christian Hof in Rscript_PresentDay.r


# Script by Susanne Fritz, last changed 13th October 2015


########################################################################

# Figure 1: the sampling map for NPP data and mammals in grid cells


# needs "Data_RegionsGridded.shp"
# needs "Data_GenusOccurrence.txt"
# needs "Data_Paleoclimate_Regions.txt"

# first, make a gridded dataset of raw diversity for each global time bin
# this only got used to plot Fig. 1.

setwd("D://Data//Gammadiv//FilesForUpload")
library(sp)
library(maps)
library(maptools)
# load the gridded regions as a spatial object
mygrid <- readShapePoly("Data_RegionsGridded.shp")
head(mygrid@data)
# load the locality-grid data
ds <- read.delim("Data_GenusOccurrence.txt", as.is=TRUE)
head(ds)
# Count the number of genera in each grid cell for each global time bin
luni <- function(x) (length(unique(x))) # a function to count unique genera in each cell or time bin
ds$WM_ID <- as.factor(ds$WM_ID) # need to make the grid cell IDs a factor to keep all of them in the tapply() bit
ds$latlong <- paste(ds$LAT, ds$LONG, sep="_") # unique location identifier
# loop through the global time bins and make a dataset where each line is a grid cell
unitime <- unique(ds$GlobalTimeBins)
for (i in seq(along=unitime)) {
  dsi <- ds[which(ds$GlobalTimeBins==unitime[i]),] # the dataset for time bin i
  gennums <- tapply(dsi$Genus, dsi$WM_ID, luni) # get the genus numbers within each grid cell
  locnums <- tapply(dsi$latlong, dsi$WM_ID, luni) # get the number of locations within each grid cell
  # collate the data into a new dataframe
  griddsi <- data.frame(WM_ID=names(gennums), GlobalTimeBins=unitime[i], NumLocs=locnums, NumGenera=gennums)
  if (i==1) (gridds <- griddsi) else (gridds <- rbind(gridds, griddsi))
}
# append the centroid information for the grid cells
gridds$Centroid_X <- ds$Centroid_X[match(gridds$WM_ID, ds$WM_ID)]
gridds$Centroid_Y <- ds$Centroid_Y[match(gridds$WM_ID, ds$WM_ID)]
# check the dataset
head(gridds)
summary(as.factor(gridds$GlobalTimeBins))
table(gridds$NumLocs, gridds$GlobalTimeBins)
length(unique(gridds$WM_ID)) # 429 grid cells with some mammal localities
tapply(gridds$WM_ID, gridds$GlobalTimeBins, luni)
summary(mygrid)
length(unique(mygrid$WM_ID)) # 1635 are in the overall grid
# make a dataset for the grid cells without mammal localities
missg <- mygrid@data$WM_ID[which(!is.element(mygrid$WM_ID, gridds$WM_ID))]
missgd <- data.frame(WM_ID=rep(missg, 8), GlobalTimeBins=rep(unitime, each=length(missg)))
summary(as.factor(missgd$GlobalTimeBins))
missgd$NumLocs <- as.numeric(NA)
missgd$NumGenera <- as.numeric(NA)
missgd$Centroid_X <- mygrid$Centroid_X[match(missgd$WM_ID, mygrid$WM_ID)]
missgd$Centroid_Y <- mygrid$Centroid_Y[match(missgd$WM_ID, mygrid$WM_ID)]
names(missgd)
names(gridds)
gridds$WM_ID <- as.numeric(as.character(gridds$WM_ID))
# append to the grid dataset and save
gridds2 <- rbind(gridds, missgd)
str(gridds2)
summary(as.factor(gridds2$GlobalTimeBins))
length(unique(gridds2$WM_ID))
write.table(gridds2, "Data_GriddedDiversity.txt", sep="\t", row.names=FALSE, col.names=TRUE)


# then, actually plot the figure

library(sp)
library(maps)
library(maptools)
# load the gridded regions as a spatial object
mygrid <- readShapePoly("Data_RegionsGridded.shp")
head(mygrid@data)
# aggregate the grid cells within each region to get region outlines
table(as.factor(mygrid$Region))
myregions <- unionSpatialPolygons(mygrid, as.factor(mygrid$Region))
plot(mygrid, border="grey")
map(add=TRUE)
plot(myregions, border="red", add=TRUE)
# load the gridded data on raw diversity
ds <- read.delim("Data_GriddedDiversity.txt", as.is=TRUE)
head(ds)
summary(as.factor(ds$GlobalTimeBins))
length(unique(ds$WM_ID))

# get the number of locations per grid cells across all time intervals
ds$NumLocs[which(is.na(ds$NumLocs))] <- 0 # turn NA into 0 so that sum works later
head(ds)
# get the number of locations for each grid cell as sum of locations within time bins
numlocs <- tapply(ds$NumLocs, ds$WM_ID, sum)
str(numlocs)
# make a new data frame with one row for each grid cell
newds <- data.frame(WM_ID=names(numlocs), NumLocs=numlocs)
newds$Centroid_X <- ds$Centroid_X[match(newds$WM_ID, ds$WM_ID)]
newds$Centroid_Y <- ds$Centroid_Y[match(newds$WM_ID, ds$WM_ID)]
length(unique(newds$WM_ID))
length(newds$WM_ID)
coordinates(newds) <- c("Centroid_X", "Centroid_Y") # turn this into a spatial dataset
gridded(newds) <- TRUE # ...and into a raster grid thingy
# make different colours from the number of locations for the colour ramping
hist(newds$NumLocs)
# log of number of locations for trying different gray ramp
newds$lognum <- log(newds$NumLocs)
newds$lognum[which(newds$NumLocs==0)] <- NA
hist(newds$lognum)
# simple number of locations between 0 and 1
newds$NumLocsProp <- 1 - (newds$NumLocs - min(newds$NumLocs)) / (max(newds$NumLocs) - min(newds$NumLocs))
hist(newds$NumLocsProp)
# ramp the logged number of locations between 0 and 1
newds$LogProp <- 1 - (newds$lognum - min(newds$lognum, na.rm=TRUE)) / (max(newds$lognum, na.rm=TRUE)-min(newds$lognum, na.rm=TRUE))
hist(newds$LogProp)
summary(newds)

# get the point locations for the NPP data

clim <- read.delim("Data_Paleoclimate_Regions.txt", as.is=TRUE)
head(clim)
summary(as.factor(clim$Source))
summary(clim$WM_ID)
summary(as.factor(clim$Continent))
summary(as.factor(clim$Region))
length(unique(clim$Spatial_ID))
luni <- function(x) (length(unique(x)))
tapply(clim$Spatial_ID, clim$Continent, luni)
summary(as.factor(clim$BestStage)) # 439 data points for best dating
summary(as.factor(clim$AutomaticStage)) # 423 data points for automatic dating

# how many locations have both MAT and MAP and were used in the calculation?
# to find this out, have to use grid cell information: if a grid cell contained any
# MAT and any MAP records within the same dating scheme, the location was used.
wmbest <- unique(clim$WM_ID[which(!is.na(clim$BestStage))]) # 100 grid cells with best-dates information
wmaut <- unique(clim$WM_ID[which(!is.na(clim$AutomaticStage))]) # 93 grid cells with automatic-dates information
all(is.element(wmaut, wmbest))
# all grid cells with automatic-dates info also have information in the
# best-dates dataset. Hence use this from now.
clim <- clim[which(!is.na(clim$BestStage)),]
# now loop through the grid cells and extract whether both MAT and MAP
# exist for each cell, and record it in the dataset.
clim$Used <- "no"
for (i in seq(along=wmbest)) {
  climi <- clim[which(clim$WM_ID==wmbest[i]),] # dataset for grid cell i
  # # which time bins are covered in any data in grid cell i?
  binsi <- unique(climi$BestStage)
  for (j in seq(along=binsi)) { # loop through these bins
    climij <- climi[which(climi$BestStage==binsi[j]),]
    # mark the locations which are in a grid cell and time bin that has data for 
    # both MAT and MAP; and chuck out the locations without both information bits
    if (any(!is.na(climij$MAT_C)) & any(!is.na(climij$MAP_mm_yr))) {
      clim$Used[which(clim$WM_ID==wmbest[i] & clim$BestStage==binsi[j])] <- "YES"
    } else {
      clim$Used[which(clim$WM_ID==wmbest[i] & clim$BestStage==binsi[j])] <- "NO"
    }
  }
}
summary(as.factor(clim$Used)) # 391 data-by-location combinations
clim <- clim[which(clim$Used=="YES"),]
summary(as.factor(clim$Continent))
summary(as.factor(clim$BestStage))
length(unique(clim$Data_ID)) # 306 unique data points
length(unique(clim$Spatial_ID)) # 145 unique spatial locations
tapply(clim$Spatial_ID, clim$Continent, luni)

# make a localities dataset
spds <- clim[which(!duplicated(clim$Spatial_ID)),]
length(spds[,1]) 

# plot the maps

pdf("Fig1.pdf", width=8, height=5)
par(mar=c(0,0,0,0))
# North America
map(col="white", fill=TRUE, lty=0, bg="lightcyan1", xlim=c(-128,-68), ylim=c(20,58))
plot(mygrid, border="light grey", add=TRUE)
image(newds[,"LogProp"], zlim=c(0,1), col = grey(seq(0.3,0.8,by=0.01)), add=TRUE)
map(add=TRUE, col="dark grey")
plot(myregions, border="black", add=TRUE, lwd=1)
points(latitude ~ longitude, spds, pch=23, col="black", bg="forestgreen", cex=0.8)
box(col="dark grey")
# Europe
map(col="white", fill=TRUE, lty=0, bg="lightcyan1", xlim=c(-10,51), ylim=c(25,58))
plot(mygrid, border="light grey", add=TRUE)
image(newds[,"LogProp"], zlim=c(0,1), col = grey(seq(0.3,0.8,by=0.01)), add=TRUE)
map(add=TRUE, col="dark grey")
plot(myregions, border="black", add=TRUE, lwd=1)
points(latitude ~ longitude, spds, pch=23, col="black", bg="forestgreen", cex=0.8)
box(col="dark grey")
# plot the legend
par(mfrow=c(1,2), mar=c(0,0,0,5))
plot(0:1, 0:1, type="n", xlab="", ylab="", xaxt="n", yaxt="n", bty="n")
# ramp of grey as in the image plots above
points(x=rep(1, 1000), y=seq(0.8,0.3,length.out=1000), pch="-", col=grey(seq(0.3,0.8,length.out=1000)))
# label with original number of localities
axis(4, at=(seq(0, 1, by=0.2)*0.5)+0.3, labels=round(exp(seq(0, 1, by=0.2)*(max(newds$lognum, na.rm=TRUE)-min(newds$lognum, na.rm=TRUE))+min(newds$lognum, na.rm=TRUE)),0), las=1)
mtext(side=4, at=0.55, text="Number of localities", line=2.5)
dev.off()


##################################################################


# Figure 2: plotted in 2 parts
# - temporal patterns of gamma diversity within continents and regions
# - temporal patterns of NPP estimates within continents

# first part:
# needs Data_Continents and Data_Regions
# needs Data_AnalysesContinents and Data_AnalysesRegions (for present-day data)

setwd("D://Data//Gammadiv//FilesForUpload")

# continental data - fossils
genall <- read.delim("Data_Continents.txt", as.is=TRUE)
head(genall)
unique(genall$Continent)
genall[which(genall$NumLocs<6),] # only MN01 in Europe has too few locations
# sort the dataset
genall <- genall[rev(order(genall$Continent, genall$GlobalBinStart, genall$ContinBinStart)),]
# make a vector with short names for the global time bins (stratigraphic stages)
summary(as.factor(genall$GlobalTimeBins))
genall$ShortBin <- substr(genall$GlobalTimeBins, 1, 2) # take the first two letters
genall$ShortBin[genall$ShortBin=="Ea"] <- "EP"
genall$ShortBin[genall$GlobalTimeBins=="Late Pliocene"] <- "LP"
summary(as.factor(genall$ShortBin))
myshortbins <- c(unique(genall$ShortBin)[-1], "Pl") # make the vector and sort it by temporal order, plus get rid of NA and add Pleistocene
# make a similar vector for the end times of bins,
# for labelling the stages into boxes along the x axis
mybinends <- c(unique(genall$GlobalBinEnd[!is.na(genall$GlobalBinEnd)]), 0.01)

# continental present-day data
dsp <- read.delim("Data_AnalysesContinents.txt", as.is=TRUE)
head(dsp)
dsp <- dsp[which(dsp$GlobalTimeBins=="Present"),]

# regional data - fossils
genallr <- read.delim("Data_Regions.txt", as.is=TRUE)
head(genallr)
unique(genallr$Region)
genallr[which(genallr$NumLocs<6),1:3] # quite a lot of region-stage combinations with not enough locations, especially land mammal ages (continental time bins)!
# sort the dataset
genallr <- genallr[rev(order(genallr$Continent, genallr$Region, genallr$GlobalBinStart, genallr$ContinBinStart)),]
# regional data - present-day
dspr <- read.delim("Data_AnalysesRegions.txt", as.is=TRUE)
head(dspr)
dspr <- dspr[which(dspr$GlobalTimeBins=="Present"),]


# a plot function specially made for the dataset
# this only works by adding to an existing plot! Also the column names need to be
# specified in quotes. The function loops through the mytime time interval names and
# matches them up to data in ds, leading to a stepped temporal trajectory to be
# plotted with error bars.
# ds - the dataset holding all columns below
# mytime - a vector of unique time interval names matching the entries in var.time
# var.time - a column name of ds, for the column that hold the time interval names
# var.time.start - a column name of ds, for the column that hold the age of the time interval start
# var.time.centre - a column name of ds, for the column that hold the age of the time interval centre
# var.time.end - a column name of ds, for the column that hold the age of the time interval end
# var.y - a column name of ds, the variable to be plotted as a stepped line against time
# var.stderr - a column name of ds, the data for error bars (e.g. std error values for the values given in var.y)
# var.sampling - a column name of ds, holding the number of sampling units; specified so data under a certain sampling level can be omitted from the plot
# sampling.min - the minimum value given in var.sampling for var.y to be plotted for a particular time interval in var.time
# mycol - the colour for the stepped line and error bars
# mylwd - the line width for the stepped line and error bars

temp.traj <- function(ds, mytime, var.time, var.time.start, var.time.centre, var.time.end, var.y, var.stderr, var.sampling, sampling.min=6, mycol="black", mylwd=1) {
	par(lend="square", lwd=mylwd) # setting some plotting parameters
	for (i in seq(along=mytime)) { # loop through the given time intervals
	    dsi <- ds[which(ds[,var.time]==mytime[i]),] # the dataset for interval i
	    ynow <- dsi[,var.y] # the value of var.y at time i
	    # if the sampling is above the sampling.min threshold and
	    # there is a var.y value at time i:
	    if (dsi[,var.sampling]>=sampling.min & !is.na(ynow)) {
	    	# plot the line for time i
	    	lines(x=c(dsi[,var.time.start], dsi[,var.time.end]), y=rep(ynow, 2), col=mycol)
	    	# plot the error bar for time i
			arrows(x0=dsi[,var.time.centre], y0=ynow+dsi[,var.stderr], x1=dsi[,var.time.centre], y1=ynow-dsi[,var.stderr], length=0, code=3, col=mycol)
			# if time i is not the first interval
			if (i > 1) {
				# get the var.y value of the previous time interval
        		ybefore <- ds[which(ds[var.time]==mytime[i-1]),var.y]
        		# if that value was not NA and sampling was good enough, plot a 
        		# vertical line to link the value.y of the previous time interval
        		# with that of time i
        		if (ds[which(ds[,var.time]==mytime[i-1]),var.sampling]>=sampling.min & !is.na(ybefore)) (lines(x=rep(dsi[var.time.start], 2), y=c(ybefore, ynow), col=mycol))
      		} # if time i is the first interval, no vertical line is needed
    	} # if there is no var.y value at time i, no lines are needed for it
  	} # close the loop through time intervals
} # close the function


# set parameters needed
my.c <- unique(genall$Continent) # list of continents
my.r <- unique(genallr$Region)[c(1,3,2,4:6)] # list of regions
# set up axis extents for diversity
summary(genall$jack1[which(genall$NumLocs>5)])
max(genall$jack1[which(genall$NumLocs>5)]+genall$jack1.se[which(genall$NumLocs>5)], na.rm=TRUE)
summary(genallr$jack1[which(genallr$NumLocs>5)])
max(genallr$jack1[which(genallr$NumLocs>5)]+genallr$jack1.se[which(genallr$NumLocs>5)], na.rm=TRUE)
dsp$jack1
dspr$jack1
miny <- 0 # y axis start for all plots
maxy <- 300 # y axis end for continental plots
maxyr <- 210 # y axis end for regional plots


# plot the figure

pdf("Fig2A-H.pdf", width=7, height=7)
par(mfcol=c(4,2), oma=c(4,3.2,1,0), mar=c(0.5,0.5,0.5,0.5), lend="square", lwd=0.9)

# loop through the continents
for (j in seq(along=my.c)) { 
  
	# continental data
	# set up plot with grey bars for every second global time bin
	plot(jack1 ~ GlobalBinCenter, genall, xlim=c(24,-0.5), ylim=c(miny, maxy), type="n", xlab="", ylab="", xaxt="n", yaxt="n", xaxs="i", yaxs="i", las=1, bty="n")
	greybins <- c("Aquitanian", "Langhian", "Tortonian", "Early Pliocene")
	for (i in seq(along=greybins)) {
		binstart <- unique(genall$GlobalBinStart[which(!is.na(genall$GlobalBinStart) & genall$GlobalTimeBins==greybins[i])])
		binend <- unique(genall$GlobalBinEnd[which(!is.na(genall$GlobalBinEnd) & genall$GlobalTimeBins==greybins[i])])
		polygon(x=c(binstart, binend, binend, binstart), y=c(rep(0,2),rep(maxy, 2)), col=grey(0.92), border=FALSE)
	}
	polygon(x=c(1.8, 0.1, 0.1, 1.8), y=c(rep(0,2),rep(maxy, 2)), col=grey(0.92), border=FALSE) # add Pleistocene
	
	# extract the data for continent j
	genallj <- genall[which(genall$Continent==my.c[j]),]
	my.r.j <- my.r[is.element(my.r, genallr$Region[which(genallr$Continent==my.c[j])])]
	
	# get the time bin info and plot diversity in continental time bins
	if (length(grep("America", my.c[j]))>0) (unicont <- unique(genall$ContinentalBins[genall$Continent==my.c[j] & !is.na(genall$ContinentalBins)])) else (unicont <- unique(genall$ContinentalBins[genall$Continent=="Europe" & !is.na(genall$ContinentalBins)]))
	temp.traj(genallj, mytime=unicont, var.time="ContinentalBins", var.time.start="ContinBinStart", var.time.centre="ContinBinCenter", var.time.end="ContinBinEnd", var.y="jack1", var.stderr="jack1.se", var.sampling="NumLocs", sampling.min=6, mycol="red3", mylwd=1.2)
  
	# diversity in the global time intervals as a smoothed line
	globalds <- genallj[which(is.na(genallj$ContinentalBins)),]
	lines(spline(globalds$GlobalBinCenter, y=globalds$jack1, method="natural", n=50), col="black", lwd=1.2)
	arrows(x0=globalds$GlobalBinCenter, y0=globalds$jack1 + globalds$jack1.se, x1=globalds$GlobalBinCenter, y1=globalds$jack1-globalds$jack1.se, length=0, code=3, lwd=1.2, col="black")
	# add the present
	points(0, dsp$jack1[dsp$Continent==my.c[j]], pch=15, col="lightskyblue", cex=1.2)
	
	# labelling and legend
	box()
	mtext(side=3, my.c[j], line=0.5, cex=0.7) # label the continent at the top
	if (j==1) { # add full y axis only for the first plot
    	# legend(0, 290, legend=c("global stratigraphic stages", "continental land mammal ages", "present day"), xjust=1, text.col=c("black", "red3", "lightskyblue"), bty="n")
    	axis(2, las=1, hadj=0.8) # y axis with labels for the left column of figures
    	text(24.2, 275, "a", pos=4, font=2, cex=1.2) # label multipart figure
	} else {
		axis(2, labels=FALSE) # y axis without labels
		text(24.2, 275, "e", pos=4, font=2, cex=1.2) # label multipart figure
	}
	axis(1, at=seq(25,0,by=-5), labels=FALSE) # x axis without labels

	# regional data
	for (k in seq(along=my.r.j)) { # loop through the regions in continent j
		# get the dataset for region k in continent j
		genallrk <- genallr[which(genallr$Region==my.r.j[k]),]
		
		# set up plot with grey bars for every second global time bin
		plot(jack1 ~ GlobalBinCenter, genallr, xlim=c(24,-0.5), ylim=c(miny, maxyr), type="n", xlab="", ylab="", xaxt="n", yaxt="n", xaxs="i", yaxs="i", bty="n")
		for (i in seq(along=greybins)) {
			binstart <- unique(genallr$GlobalBinStart[which(!is.na(genallr$GlobalBinStart) & genallr$GlobalTimeBins==greybins[i])])
			binend <- unique(genallr$GlobalBinEnd[which(!is.na(genallr$GlobalBinEnd) & genallr$GlobalTimeBins==greybins[i])])
			polygon(x=c(binstart, binend, binend, binstart), y=c(rep(0,2),rep(maxy, 2)), col=grey(0.92), border=FALSE)
		}
		polygon(x=c(1.8, 0.1, 0.1, 1.8), y=c(rep(0,2),rep(maxy, 2)), col=grey(0.92), border=FALSE) # add Pleistocene
		
		# plot the continental diversity through time
		temp.traj(genallrk, mytime=unicont, var.time="ContinentalBins", var.time.start="ContinBinStart", var.time.centre="ContinBinCenter", var.time.end="ContinBinEnd", var.y="jack1", var.stderr="jack1.se", var.sampling="NumLocs", sampling.min=6, mycol="red3", mylwd=1.2)
		
		# diversity in the global time intervals as a smoothed line
		globalds <- genallrk[which(is.na(genallrk$ContinentalBins)),]
		lines(spline(globalds$GlobalBinCenter, y=globalds$jack1, method="natural", n=50), col="black", lwd=1.2)
		arrows(x0=globalds$GlobalBinCenter, y0=globalds$jack1 + globalds$jack1.se, x1=globalds$GlobalBinCenter, y1=globalds$jack1-globalds$jack1.se, length=0, code=3, lwd=1.2, col="black")
		# add the present
		points(0, dspr$jack1[dspr$Region==my.r.j[k]], pch=15, col="lightskyblue", cex=1.2)
		
		# labelling and legend
		if (j==1) {
			axis(2, las=1, hadj=0.8) # y axis with labels for the left column of figures
			text(24.2, 190, letters[2:4][k], pos=4, font=2, cex=1.2) # label multipart figure
		} else { 
			axis(2, labels=FALSE) # y axis without labels for the right column of figures
			text(24.2, 190, letters[6:8][k], pos=4, font=2, cex=1.2) # label multipart figure
		}
		if (k==3) { # for the bottom row of figures, make a nice x axis
			# add boxes and labels showing the stratigraphic stages 
			for (i in seq(along=c(23, mybinends))) {
				if (is.element(i, c(2,4,6,8,10))) (polygon(x=c(c(23, mybinends)[i-1], c(23, mybinends)[i], c(23, mybinends)[i], c(23, mybinends)[i-1]), y=c(rep(-20,2),rep(miny, 2)), col=grey(0.92), border=FALSE, xpd=NA))
			}
			for (i in seq(along=c(23, mybinends))) {
				lines(x=c(c(23, mybinends)[i], c(23, mybinends)[i]), y=c(-20, miny), xpd=NA)
			}
			lines(x=c(24,24), y=c(-20, miny), xpd=NA)
			lines(x=c(-0.5,-0.5), y=c(-20, miny), xpd=NA)
			lines(x=c(24,-0.5), y=c(-20,-20), xpd=NA)
			text(x=c(rev(sort(unique(genall$GlobalBinCenter[!is.na(genall$GlobalBinCenter)]))), 0.9), y=-10, myshortbins, xpd=NA)
			axis(1, pos=-20, padj=-0.5)
			mtext(side=1, "Time (million years ago)", line=3, cex=0.7)
		} else (axis(1, at=seq(25,0,by=-5), labels=FALSE))
		box()
	} # loop through the regions
	if (j==1) (mtext(expression(paste("Mammalian ", gamma, " diversity")), 2, outer=TRUE, line=2, cex=0.7)) # label the y axis for all panels
} # loop through the continents
dev.off()


# second part:
# only needs Data_Continents and Data_AnalysesContinents (for present-day data)

# how many entries are based on one grid cell only?
genall[which(genall$NPP_AutoPlants_NumCells<2),]
genall[which(genall$NPP_BestPlants_NumCells<2),]
# size of NPP points shows number of data points (ie grid cells with data)
genall$SymbolSizeAut <- 0.5
genall$SymbolSizeAut[which(genall$NPP_AutoPlants_NumCells>2 & genall$NPP_AutoPlants_NumCells<6)] <- 0.75
genall$SymbolSizeAut[which(genall$NPP_AutoPlants_NumCells>5 & genall$NPP_AutoPlants_NumCells<16)] <- 1
genall$SymbolSizeAut[which(genall$NPP_AutoPlants_NumCells>15)] <- 1.3
genall$SymbolSizeAut[which(genall$NPP_AutoPlants_NumCells<2)] <- NA
genall$SymbolSizeBest <- 0.5
genall$SymbolSizeBest[which(genall$NPP_BestPlants_NumCells>2 & genall$NPP_BestPlants_NumCells<6)] <- 0.75
genall$SymbolSizeBest[which(genall$NPP_BestPlants_NumCells>5 & genall$NPP_BestPlants_NumCells<16)] <- 1
genall$SymbolSizeBest[which(genall$NPP_BestPlants_NumCells>15)] <- 1.3
genall$SymbolSizeBest[which(genall$NPP_BestPlants_NumCells<2)] <- NA
# set parameters needed
my.c <- c("North America", "Europe")
# limits of the NPP data (and axis extents for the plot)
miny <- min(c(genall$NPP_BestPlants_meanminwght50, genall$NPP_AutoPlants_meanminwght50, dsp$NPP_meanhdrwght50-dsp$NPP_hdrstderr), na.rm=TRUE) - 120
maxy <- max(c(genall$NPP_BestPlants_meanmaxwght50, genall$NPP_AutoPlants_meanmaxwght50, dsp$NPP_meanhdrwght50+dsp$NPP_hdrstderr), na.rm=TRUE) + 40


pdf("Fig2J-K.pdf", width=7, height=2.5)
par(mfrow=c(1,2), lend="square", oma=c(2.7,5,1,0), mar=c(0.5,0.5,0.5,0.5), cex=0.7)

# loop through the continents
for (j in seq(along=my.c)) { 
  climc <- genall[which(genall$Continent==my.c[j]),] # select data for continent j
  
  # set up plot with grey bars for every second global time bin
  plot(NPP_BestPlants_meanhdrwght50 ~ GlobalBinCenter, genall, xlim=c(24,-0.5), ylim=c(miny, maxy), type="n", xlab="", ylab="", xaxt="n", yaxt="n", xaxs="i", yaxs="i", las=1, bty="n")
  greybins <- c("Aquitanian", "Langhian", "Tortonian", "Early Pliocene")
  for (i in seq(along=greybins)) {
    binstart <- unique(genall$GlobalBinStart[which(!is.na(genall$GlobalBinStart) & genall$GlobalTimeBins==greybins[i])])
    binend <- unique(genall$GlobalBinEnd[which(!is.na(genall$GlobalBinEnd) & genall$GlobalTimeBins==greybins[i])])
    polygon(x=c(binstart, binend, binend, binstart), y=c(rep(0,2),rep(maxy, 2)), col=grey(0.92), border=FALSE)
  }
  polygon(x=c(1.8, 0.1, 0.1, 1.8), y=c(rep(0,2),rep(maxy, 2)), col=grey(0.92), border=FALSE) # add Pleistocene
  
  # plot NPP data - best dates
  lines(climc$GlobalBinCenter+0.25, climc$NPP_BestPlants_meanhdrwght50, col="forestgreen", lty=2) # background line to easily follow temporal trajectory
  arrows(x0=climc$GlobalBinCenter+0.25, y0=climc$NPP_BestPlants_meanminwght50, x1=climc$GlobalBinCenter+0.25, y1=climc$NPP_BestPlants_meanmaxwght50, length=0, code=3, col="forestgreen") # error bars
  points(climc$GlobalBinCenter+0.25, climc$NPP_BestPlants_meanhdrwght50, pch=16, col="forestgreen", cex=climc$SymbolSizeBest) # means for stages with 2 or more grid cells
  if (any(climc$NPP_BestPlants_NumCells[!is.na(climc$NPP_BestPlants_NumCells)]==1)) (points(climc$GlobalBinCenter[climc$NPP_BestPlants_NumCells==1]+0.25, climc$NPP_BestPlants_meanhdrwght50[climc$NPP_BestPlants_NumCells==1], pch=4, cex=0.7, col="forestgreen")) # means for stages with one grid cell
  
  # plot NPP data - automatic dates
  arrows(x0=climc$GlobalBinCenter-0.25, y0=climc$NPP_AutoPlants_meanminwght50, x1=climc$GlobalBinCenter-0.25, y1=climc$NPP_AutoPlants_meanmaxwght50, length=0, code=3, col="orange") # error bars
  points(climc$GlobalBinCenter-0.25, climc$NPP_AutoPlants_meanhdrwght50, pch=16, col="orange", cex=climc$SymbolSizeAut) # means for stages with 2 or more grid cells
  if (any(climc$NPP_AutoPlants_NumCells[!is.na(climc$NPP_AutoPlants_NumCells)]==1)) (points(climc$GlobalBinCenter[climc$NPP_AutoPlants_NumCells==1]-0.25, climc$NPP_AutoPlants_meanhdrwght50[climc$NPP_AutoPlants_NumCells==1], pch=4, cex=0.7, col="orange")) # means for stages with one grid cell
  
  # add the present points with error bars
  dspc <- dsp[which(dsp$Continent==my.c[j]),] # select present data for continent j
  arrows(x0=dspc$GlobalBinCenter, y0=dspc$NPP_meanhdrwght50 - dspc$NPP_hdrstderr, x1=dspc$GlobalBinCenter, y1=dspc$NPP_meanhdrwght50 + dspc$NPP_hdrstderr, length=0, code=3, col="lightskyblue") # error bar
  points(NPP_meanhdrwght50 ~ GlobalBinCenter, dspc, pch=15, col="lightskyblue") # mean estimate
  
  # labelling and legend
  box()
  # label the continents at the top
  mtext(my.c[j], side=3, line=0.2, cex=0.8)
  # sort out the axes and axis labels
  if (j==1) {
    axis(2, las=1, cex.axis=0.9)
    mtext("Net primary production", 2, line=4.3, cex=0.8)
    mtext(expression(paste("(g dry matter/m"^2, "/year)")), 2, line=2.8, cex=0.8)
  } else (axis(2, labels=FALSE))
  axis(1, at=seq(25,0,by=-5), padj=-0.5, cex.axis=0.9) # x axis
  # add boxes and labels showing the stratigraphic stages 
  for (i in seq(along=c(23, mybinends))) {
    lines(x=c(c(23, mybinends)[i], c(23, mybinends)[i]), y=c(miny, miny+80))
  }
  lines(x=c(24,-0.5), y=c(miny+80, miny+80))
  text(x=c(unique(genall$GlobalBinCenter[!is.na(genall$GlobalBinCenter)]), 0.9), y=miny+40, myshortbins)
  mtext(side=1, "Time (million years ago)", line=1.9, cex=0.8)
  # legend for sample size (number of grid cells used)
  if (j==1) {
    legend(22.8, 1800, title="sample size", legend=c("1", "2","3-5","6-15","16-34"), col="black", pch=c(4, rep(16, 4)), cex=0.9, pt.cex=c(0.7,0.5,0.75,1,1.3), bty="n")
  }
}
dev.off()


##########################################################

# Figure 3: diversity-climate models and relationships
# plus comparison with the present day

# needs Data_AnalysesContinents and Data_AnalysesRegions

setwd("D://Data//Gammadiv//FilesForUpload")
library(lme4)

# load all datasets in: Continents
dats <- read.delim("Data_AnalysesContinents.txt", as.is=TRUE)
dats$Continent <- as.factor(dats$Continent)
str(dats)
datsc <- dats
table(datsc$Continent, datsc$GlobalTimeBins)
# Regions
dats <- read.delim("Data_AnalysesRegions.txt", as.is=TRUE)
dats$Region <- as.factor(dats$Region)
str(dats)
datsr <- dats
table(datsr$Region, datsr$GlobalTimeBins)

pdf("Fig3.pdf", width=8, height=4)
par(mfrow=c(1,2), mar=c(5,1,0.5,0.5), oma=c(0,3,0,0))

# CONTINENTS
dats <- datsc
dats$ShortBin <- substr(dats$GlobalTimeBins, 1, 2)
dats$ShortBin[dats$ShortBin=="Ea"] <- "EP"
# for the y axis, get all response values to make the same y axis
# for continents and regions, plus add 25 as the minimum value for the present
myy <- c(dats$BestModelResiduals, datsr$BestModelResiduals, dats$jack1[dats$GlobalTimeBins=="Present"], dats$jack1_adj_extinction[dats$GlobalTimeBins=="Present"], datsr$jack1[datsr$GlobalTimeBins=="Present"], datsr$jack1_adj_extinction[datsr$GlobalTimeBins=="Present"])
# the NPP values for the x axis also for both continents and regions, with 470 for the present day
myx <- c(dats$NPP_meanhdrwght50, datsr$NPP_meanhdrwght50, dats$NPP_adj_HANPP[dats$GlobalTimeBins=="Present"], datsr$NPP_adj_HANPP[dats$GlobalTimeBins=="Present"])
# set up continental plot
plot(BestModelResiduals ~ NPP_meanhdrwght50, dats, type="n", xaxt="n", yaxt="n", las=1, xlim=c(min(myx, na.rm=TRUE), max(myx, na.rm=TRUE)), ylim=c(min(myy, na.rm=TRUE)-5, max(myy, na.rm=TRUE)+12), xlab="", ylab="")
axis(1)
mtext("Net primary production", 1, line=2.3, cex=1)
mtext(expression(paste("(g dry matter/m"^2, "/year)")), 1, line=3.4, cex=1)
axis(2, las=1)
mtext(expression(paste("Mammalian ", gamma, " diversity")), 2, line=2.8, cex=1)

# plot the model line: get the right model
dats$jack1.int <- round(dats$jack1, 0)
mymodel <- glmer(jack1.int ~ NPP_meanhdrwght50_lg + Area_lg + Time_lg + (GlobalBinCenter | Continent), dats[which(dats$GlobalTimeBins!="Present"),], family=poisson)
summary(mymodel)
# make a dataframe for the fixed effects
newdat <- expand.grid(Myvar=seq(log(min(myx, na.rm=TRUE)), log(max(myx, na.rm=TRUE)), length.out=2000), Area=median(dats$Area_lg, na.rm=TRUE), Time=median(dats$Time_lg, na.rm=TRUE), Resp=0)
names(newdat) <- c("NPP_meanhdrwght50_lg", "Area_lg", "Time_lg", "jack1.int")
# predict the response variable for the fixed effects
mm <- model.matrix(terms(mymodel), newdat)
newdat[,length(newdat)] <- exp(mm %*% getME(mymodel, "fixef"))
# plot the line: continuous for the fossil data and predict outside range as dashed line
minfossil <- min(dats$NPP_meanhdrwght50_lg[dats$GlobalTimeBins!="Present"], na.rm=TRUE) # the minimum NPP for just the fossil data points
lines(exp(newdat[which(newdat$NPP_meanhdrwght50_lg >= minfossil),"NPP_meanhdrwght50_lg"]), newdat[which(newdat$NPP_meanhdrwght50_lg >= minfossil),length(newdat)], col="black", lwd=1.5)
lines(exp(newdat[which(newdat$NPP_meanhdrwght50_lg < minfossil),"NPP_meanhdrwght50_lg"]), newdat[which(newdat$NPP_meanhdrwght50_lg < minfossil),length(newdat)], col="black", lty=2, lwd=1.5)
# plot an upper and lower line to show random effects
fixr <- getME(mymodel, "fixef") # original parameter estimates
vcm <- as.data.frame(VarCorr(mymodel)) # variance estimates of random effects
fixr[1] <- fixr[1] + vcm[1,4] # add variance to intercept
fixr[2] <- fixr[2] + vcm[2,4] # add variance to slope
newdat[,length(newdat)] <- exp(mm %*% fixr) # predict the upper estimated relationship for random effect variance
lines(exp(newdat$NPP_meanhdrwght50_lg[which(newdat$NPP_meanhdrwght50_lg >= minfossil)]), newdat[which(newdat$NPP_meanhdrwght50_lg >= minfossil),length(newdat)], lty=3, lwd=1.5)
fixr[1] <- fixr[1] - 2*vcm[1,4] # subtract twice the variance for intercept
fixr[2] <- fixr[2] - 2*vcm[2,4] # subtract twice the variance for slope
newdat[,length(newdat)] <- exp(mm %*% fixr) # predict the lower estimated relationship for random effect variance
lines(exp(newdat$NPP_meanhdrwght50_lg[which(newdat$NPP_meanhdrwght50_lg >= minfossil)]), newdat[which(newdat$NPP_meanhdrwght50_lg >= minfossil),length(newdat)], lty=3, lwd=1.5)

# add fossil data points (conditional model responses)
points(BestModelResiduals ~ NPP_meanhdrwght50, dats[which(dats$Continent=="Europe"),], pch=16, col="black", cex=2)
points(BestModelResiduals ~ NPP_meanhdrwght50, dats[which(dats$Continent=="North America"),], pch=15, col="black", cex=2)
text(y=dats$BestModelResiduals, x=dats$NPP_meanhdrwght50, labels=dats$ShortBin, cex=0.6, col="white")
# points for extant observed diversity and HANPP and extinction adjustments
points(jack1 ~ NPP_meanhdrwght50, dats[which(dats$Continent=="Europe" & dats$GlobalTimeBins=="Present"),], pch=16, col="light grey", cex=1.8)
points(jack1 ~ NPP_meanhdrwght50, dats[which(dats$Continent=="North America" & dats$GlobalTimeBins=="Present"),], pch=15, col="light grey", cex=1.8)
points(jack1 ~ NPP_adj_HANPP, dats[which(dats$Continent=="North America" & dats$GlobalTimeBins=="Present"),], pch=15, col="lightskyblue", cex=1.8)
points(jack1 ~ NPP_adj_HANPP, dats[which(dats$Continent=="Europe" & dats$GlobalTimeBins=="Present"),], pch=16, col="lightskyblue", cex=1.8)
points(jack1_adj_extinction ~ NPP_meanhdrwght50, dats[which(dats$Continent=="Europe" & dats$GlobalTimeBins=="Present"),], pch=16, col="red3", cex=1.8)
points(jack1_adj_extinction ~ NPP_meanhdrwght50, dats[which(dats$Continent=="North America" & dats$GlobalTimeBins=="Present"),], pch=15, col="red3", cex=1.8)

# add a legend
legend(600, 265, legend=c("North America", "Europe"), pch=c(0,1), col="black", pt.cex=1, cex=0.63, bg=NA, bty="n", y.intersp=1.3, x.intersp=1.2)
text(490, 253, "a", font=2, las=1)
# text(460, 244, labels="fossil", pos=4, cex=0.7)
# text(460, 232, labels="present observed", pos=4, cex=0.7, col="grey")
# text(460, 220, labels="present adjusted for HANPP", pos=4, cex=0.7, col="lightskyblue")
# text(460, 208, labels="present adjusted for extinctions", pos=4, cex=0.7, col="red3")


# REGIONS

dats <- datsr
dats$ShortBin <- substr(dats$GlobalTimeBins, 1, 2)
dats$ShortBin[dats$ShortBin=="Ea"] <- "EP"
# set up regional plot with same axis limits as continental plot
plot(BestModelResiduals ~ NPP_meanhdrwght50, dats, type="n", xaxt="n", yaxt="n", las=1, xlim=c(min(myx, na.rm=TRUE), max(myx, na.rm=TRUE)), ylim=c(min(myy, na.rm=TRUE)-5, max(myy, na.rm=TRUE)+12), xlab="", ylab="")
axis(1)
mtext("Net primary production", 1, line=2.3, cex=1)
mtext(expression(paste("(g dry matter/m"^2, "/year)")), 1, line=3.4, cex=1)
axis(2, labels=FALSE)

# plot the model line: get the right model
dats$jack1.int <- round(dats$jack1, 0)
mymodel <- glmer(jack1.int ~ NPP_meanhdrwght50_lg + Area_lg + Time_lg + (GlobalBinCenter | Region), dats[which(dats$GlobalTimeBins!="Present"),], family=poisson)
summary(mymodel)
# make a dataframe for the fixed effects
newdat <- expand.grid(Myvar=seq(log(min(myx, na.rm=TRUE)), log(max(myx, na.rm=TRUE)), length.out=2000), Area=median(dats$Area_lg, na.rm=TRUE), Time=median(dats$Time_lg, na.rm=TRUE), Resp=0)
names(newdat) <- c("NPP_meanhdrwght50_lg", "Area_lg", "Time_lg", "jack1.int")
# predict the response variable for the fixed effects
mm <- model.matrix(terms(mymodel), newdat)
newdat[,length(newdat)] <- exp(mm %*% getME(mymodel, "fixef"))
# plot the line: continuous for the fossil data and predict outside range as dashed line
minfossil <- min(dats$NPP_meanhdrwght50_lg[dats$GlobalTimeBins!="Present"], na.rm=TRUE) # the minimum NPP for just the fossil data points
lines(exp(newdat[which(newdat$NPP_meanhdrwght50_lg >= minfossil),"NPP_meanhdrwght50_lg"]), newdat[which(newdat$NPP_meanhdrwght50_lg >= minfossil),length(newdat)], col="black", lwd=1.5)
lines(exp(newdat[which(newdat$NPP_meanhdrwght50_lg < minfossil),"NPP_meanhdrwght50_lg"]), newdat[which(newdat$NPP_meanhdrwght50_lg < minfossil),length(newdat)], col="black", lty=2, lwd=1.5)
# plot an upper and lower line to show random effects
fixr <- getME(mymodel, "fixef") # original parameter estimates
vcm <- as.data.frame(VarCorr(mymodel)) # variance estimates of random effects
fixr[1] <- fixr[1] + vcm[1,4] # add variance to intercept
fixr[2] <- fixr[2] + vcm[2,4] # add variance to slope
newdat[,length(newdat)] <- exp(mm %*% fixr) # predict the upper estimated relationship for random effect variance
lines(exp(newdat$NPP_meanhdrwght50_lg[which(newdat$NPP_meanhdrwght50_lg >= minfossil)]), newdat[which(newdat$NPP_meanhdrwght50_lg >= minfossil),length(newdat)], lty=3, lwd=1.5)
fixr[1] <- fixr[1] - 2*vcm[1,4] # subtract twice the variance for intercept
fixr[2] <- fixr[2] - 2*vcm[2,4] # subtract twice the variance for slope
newdat[,length(newdat)] <- exp(mm %*% fixr) # predict the lower estimated relationship for random effect variance
lines(exp(newdat$NPP_meanhdrwght50_lg[which(newdat$NPP_meanhdrwght50_lg >= minfossil)]), newdat[which(newdat$NPP_meanhdrwght50_lg >= minfossil),length(newdat)], lty=3, lwd=1.5)

# add fossil data points for regions (conditional model responses)
points(BestModelResiduals ~ NPP_meanhdrwght50, dats[which(dats$Region=="Eastern Europe"),], pch=16, col="black", cex=2)
points(BestModelResiduals ~ NPP_meanhdrwght50, dats[which(dats$Region=="Western Europe"),], pch=17, col="black", cex=2)
points(BestModelResiduals ~ NPP_meanhdrwght50, dats[which(dats$Region=="Western North America"),], pch=15, col="black", cex=2)
text(y=dats$BestModelResiduals, x=dats$NPP_meanhdrwght50, labels=dats$ShortBin, cex=0.6, col="white")
# points for extant observed diversity and HANPP and extinction adjustments
points(jack1 ~ NPP_meanhdrwght50, dats[which(dats$Region=="Eastern Europe" & dats$GlobalTimeBins=="Present"),], pch=17, col="light grey", cex=1.8)
points(jack1 ~ NPP_meanhdrwght50, dats[which(dats$Region=="Western Europe" & dats$GlobalTimeBins=="Present"),], pch=16, col="light grey", cex=1.8)
points(jack1 ~ NPP_meanhdrwght50, dats[which(dats$Region=="Western North America" & dats$GlobalTimeBins=="Present"),], pch=15, col="light grey", cex=1.8)
points(jack1 ~ NPP_adj_HANPP, dats[which(dats$Region=="Eastern Europe" & dats$GlobalTimeBins=="Present"),], pch=17, col="lightskyblue", cex=1.8)
points(jack1 ~ NPP_adj_HANPP, dats[which(dats$Region=="Western North America" & dats$GlobalTimeBins=="Present"),], pch=15, col="lightskyblue", cex=1.8)
points(jack1 ~ NPP_adj_HANPP, dats[which(dats$Region=="Western Europe" & dats$GlobalTimeBins=="Present"),], pch=16, col="lightskyblue", cex=1.8)
points(jack1_adj_extinction ~ NPP_meanhdrwght50, dats[which(dats$Region=="Eastern Europe" & dats$GlobalTimeBins=="Present"),], pch=17, col="red3", cex=1.8)
points(jack1_adj_extinction ~ NPP_meanhdrwght50, dats[which(dats$Region=="Western North America" & dats$GlobalTimeBins=="Present"),], pch=15, col="red3", cex=1.8)
points(jack1_adj_extinction ~ NPP_meanhdrwght50, dats[which(dats$Region=="Western Europe" & dats$GlobalTimeBins=="Present"),], pch=16, col="red3", cex=1.8)
# add a legend
legend(600, 265, legend=c("Western North America", "Western Europe", "Eastern Europe"), pch=c(0,1,2), col="black", pt.cex=1, cex=0.63, bg=NA, bty="n", y.intersp=1.3, x.intersp=1.2)
text(490, 253, "b", font=2, las=1)
dev.off()


###########################################################


# Supplementary figure 1: raw genus counts, range-through genera, and
# number of localities

# needs Data_Continents, Data_Regions; Output_RangeThruDataContinents, Output_RangeThruDataRegions

setwd("D://Data//Gammadiv//FilesForUpload")

# continental data
genall <- read.delim("Data_Continents.txt", as.is=TRUE)
head(genall)
unique(genall$Continent)
genrt <- read.delim("Output_RangeThruDataContinents.txt", as.is=TRUE)
head(genrt)
genall <- merge(genall, genrt, by=c(1:3))
head(genall)
table(genall$Continent, genall$GlobalTimeBins)
table(genall$Continent, genall$ContinentalBins)
# make a vector with short names for the global time bins (stratigraphic stages)
summary(as.factor(genall$GlobalTimeBins))
genall$ShortBin <- substr(genall$GlobalTimeBins, 1, 2) # take the first two letters
genall$ShortBin[genall$ShortBin=="Ea"] <- "EP"
genall$ShortBin[genall$GlobalTimeBins=="Late Pliocene"] <- "LP"
summary(as.factor(genall$ShortBin))
myshortbins <- unique(genall$ShortBin[rev(order(genall$GlobalBinStart))])[-1] # make the vector and sort it by temporal order, plus get rid of NA
# make a similar vector for the end times of bins,
# for labelling the stages into boxes along the x axis
mybinends <- rev(sort(unique(genall$GlobalBinEnd[!is.na(genall$GlobalBinEnd)])))
# sort the dataset
genall <- genall[rev(order(genall$Continent, genall$GlobalBinStart, genall$ContinBinStart)),]

# regional data
genallr <- read.delim("Data_Regions.txt", as.is=TRUE)
head(genallr)
unique(genallr$Region)
genrt <- read.delim("Output_RangeThruDataRegions.txt", as.is=TRUE)
head(genrt)
genrt <- genrt[,-length(genrt)] # get rid of continent column for merging
genallr <- merge(genallr, genrt, by=c(1:3), all.x=TRUE)
head(genallr)
table(genallr$Region, genallr$GlobalTimeBins)
table(genallr$Region, genallr$ContinentalBins)
genallr <- genallr[rev(order(genallr$Region, genallr$GlobalBinStart, genallr$ContinBinStart)),]
genallr$NumGenera[which(genallr$NumLocs==0)] <- 0

# a plot function specially made for the dataset, see Figure 2 above.
# changed by taking out the sampling variable and by only adding the error bar on top, not adding and subtracting.
temp.traj <- function(ds, mytime, var.time, var.time.start, var.time.centre, var.time.end, var.y, var.stderr, mycol="black", mylwd=1) {
  par(lend="square", lwd=mylwd) # setting some plotting parameters
  for (i in seq(along=mytime)) { # loop through the given time intervals
    dsi <- ds[which(ds[,var.time]==mytime[i]),] # the dataset for interval i
    ynow <- dsi[,var.y] # the value of var.y at time i
    # if there is a var.y value at time i:
    if (!is.na(ynow)) {
      # plot the line for time i
      lines(x=c(dsi[,var.time.start], dsi[,var.time.end]), y=rep(ynow, 2), col=mycol)
      # plot the error bar for time i if it has length >0
      if (dsi[,var.stderr]>0) (arrows(x0=dsi[,var.time.centre], y0=ynow+dsi[,var.stderr], x1=dsi[,var.time.centre], y1=ynow, length=0.03, code=1, angle=90, col=mycol))
      # if time i is not the first interval
      if (i > 1) {
        # get the var.y value of the previous time interval
        ybefore <- ds[which(ds[var.time]==mytime[i-1]),var.y]
        # if that value was not NA, plot a vertical line to link the value.y
        # of the previous time interval with that of time i
        if (!is.na(ybefore)) (lines(x=rep(dsi[var.time.start], 2), y=c(ybefore, ynow), col=mycol))
      } # if time i is the first interval, no vertical line is needed
    } # if there is no var.y value at time i, no lines are needed for it
  } # close the loop through time intervals
} # close the function

# set parameters needed
my.c <- unique(genall$Continent) # list of continents
my.r <- unique(genallr$Region)[c(1,5,3,2,4,6)] # list of regions
# set up axis extents for diversity
summary(genall$NumGenera)
max(genall$NumGenera+genall$NumGRangeThru1, na.rm=TRUE)
summary(genallr$NumGenera)
max(genallr$NumGenera+genallr$NumGRangeThru1, na.rm=TRUE)
miny <- 0 # y axis start for all plots
maxy <- 230 # y axis end for continental plots
maxyr <- 160 # y axis end for regional plots
summary(genall$NumLocs)
summary(genallr$NumLocs)
genall$NL_log <- log(genall$NumLocs) * 20
genallr$NL_log <- log(genallr$NumLocs) *20
summary(genall$NL_log)
summary(genallr$NL_log)
genallr$NL_log[is.infinite(genallr$NL_log)] <- 0
myy <- 150 # continental offset of the number of location y axis to the number of genera y axis
myy.r <- 120 # regional offset

pdf("FigS1.pdf", width=7, height=8)
par(mfcol=c(4,2), oma=c(4,2.8,1,2.8), mar=c(0.5,1,0.5,1), lend="square", lwd=0.9)

# loop through the continents
for (j in seq(along=my.c)) { 
  
  # continental data
  # set up plot with grey bars for every second global time bin
  plot(NumGenera ~ GlobalBinCenter, genall, xlim=c(23.8,1.8), ylim=c(miny-myy, maxy), type="n", xlab="", ylab="", xaxt="n", yaxt="n", xaxs="i", yaxs="i", las=1, bty="n")
  greybins <- c("Aquitanian", "Langhian", "Tortonian", "Early Pliocene")
  for (i in seq(along=greybins)) {
    binstart <- unique(genall$GlobalBinStart[which(!is.na(genall$GlobalBinStart) & genall$GlobalTimeBins==greybins[i])])
    binend <- unique(genall$GlobalBinEnd[which(!is.na(genall$GlobalBinEnd) & genall$GlobalTimeBins==greybins[i])])
    polygon(x=c(binstart, binend, binend, binstart), y=c(rep(miny-myy,2),rep(maxy, 2)), col=grey(0.92), border=FALSE)
  }
  abline(h=0, lwd=1)
  
  # extract the data for continent j
  genallj <- genall[which(genall$Continent==my.c[j]),]
  my.r.j <- my.r[is.element(my.r, genallr$Region[which(genallr$Continent==my.c[j])])]
  
  # get the time bin info and plot diversity in continental time bins
  if (length(grep("America", my.c[j]))>0) (unicont <- unique(genall$ContinentalBins[genall$Continent==my.c[j] & !is.na(genall$ContinentalBins)])) else (unicont <- unique(genall$ContinentalBins[genall$Continent=="Europe" & !is.na(genall$ContinentalBins)]))
  temp.traj(genallj, mytime=unicont, var.time="ContinentalBins", var.time.start="ContinBinStart", var.time.centre="ContinBinCenter", var.time.end="ContinBinEnd", var.y="NumGenera", var.stderr="NumGRangeThru1", mycol="red3", mylwd=1.2)
  
  # diversity in the global time intervals
  globalds <- genallj[which(is.na(genallj$ContinentalBins)),]
  lines(spline(globalds$GlobalBinCenter, y=globalds$NumGenera, method="natural", n=50), col="black", lwd=1.2)
  points(NumGenera ~ GlobalBinCenter, globalds, pch=16, col="black", cex=0.8)
  arrows(x0=globalds$GlobalBinCenter, y0=globalds$NumGenera + globalds$NumGRangeThru1, x1=globalds$GlobalBinCenter, y1=globalds$NumGenera, length=0.03, angle=90, code=1, lwd=1.2, col="black")
  
  # number of localities as a bar chart
  for (i in seq(along=unique(globalds$GlobalTimeBins))) { # in global time intervals
    dsi <- globalds[which(globalds$GlobalTimeBins==unique(globalds$GlobalTimeBins)[i]),]
    if (dsi$NumLocs!=0) (polygon(x=c(dsi$GlobalBinStart, dsi$GlobalBinStart, dsi$GlobalBinEnd, dsi$GlobalBinEnd), y=c(-myy, dsi$NL_log-myy, dsi$NL_log-myy, -myy), col="lightskyblue", border="navy"))
  }
  for (i in seq(along=unicont)) { # in continental bins
    dsi <- genallj[which(genallj$ContinentalBins==unicont[i]),]
    if (dsi$NumLocs!=0) (polygon(x=c(dsi$ContinBinStart, dsi$ContinBinStart, dsi$ContinBinEnd, dsi$ContinBinEnd), y=c(-myy, dsi$NL_log-myy, dsi$NL_log-myy, -myy), col=rgb(1,1,1, alpha=0.4)))
  }
  abline(h=-myy+log(5)*20, lty=3)
  
  # labelling and legend
  box()
  mtext(side=3, my.c[j], line=0.5, cex=0.7) # label the continent at the top
  if (j==1) { # add a legend only for the first plot, and label y axes
    legend(14, 230, legend=c("global stratigraphic stages", "continental land mammal ages", "number of locations"), text.col=c("black", "red3", "blue"), bty="n", cex=0.9)
    axis(2, at=pretty(genall$NumGenera), las=1, hadj=0.8, cex.axis=0.9) # y axis with labels for the left side
    axis(4, at=c(-myy, log(c(10,50,100,200,500))*20-myy), labels=FALSE) # y axis without labels for the right side
    text(24.2, 210, "A", pos=4, font=2, cex=1.2) # label multipart figure
  } else {
    axis(2, at=pretty(genall$NumGenera), labels=FALSE) # y axis without labels on left side
    axis(4, at=c(-myy, (log(c(10,50,100,200,500))*20)-myy), labels=c(0,10,50,100,200,500), las=1, hadj=0.2, cex.axis=0.9) # y axis with labels on left side
    text(24.2, 210, "E", pos=4, font=2, cex=1.2) # label multipart figure
  }
  axis(1, at=seq(25,0,by=-5), labels=FALSE) # x axis without labels
  
  # regional data
  for (k in seq(along=my.r.j)) { # loop through the regions in continent j
    # get the dataset for region k in continent j
    genallrk <- genallr[which(genallr$Region==my.r.j[k]),]
    
    # set up plot with grey bars for every second global time bin
    plot(NumGenera ~ GlobalBinCenter, genallr, xlim=c(23.8,1.8), ylim=c(miny-myy.r, maxyr), type="n", xlab="", ylab="", xaxt="n", yaxt="n", xaxs="i", yaxs="i", bty="n")
    for (i in seq(along=greybins)) {
      binstart <- unique(genallr$GlobalBinStart[which(!is.na(genallr$GlobalBinStart) & genallr$GlobalTimeBins==greybins[i])])
      binend <- unique(genallr$GlobalBinEnd[which(!is.na(genallr$GlobalBinEnd) & genallr$GlobalTimeBins==greybins[i])])
      polygon(x=c(binstart, binend, binend, binstart), y=c(rep(-myy.r,2),rep(maxy, 2)), col=grey(0.92), border=FALSE)
    }
    abline(h=0, lwd=1)
    
    # plot the continental diversity through time
    temp.traj(genallrk, mytime=unicont, var.time="ContinentalBins", var.time.start="ContinBinStart", var.time.centre="ContinBinCenter", var.time.end="ContinBinEnd", var.y="NumGenera", var.stderr="NumGRangeThru1", mycol="red3", mylwd=1.2)
    # diversity in the global time intervals
    globalds <- genallrk[which(is.na(genallrk$ContinentalBins)),]
    lines(spline(globalds$GlobalBinCenter, y=globalds$NumGenera, method="natural", n=50), col="black", lwd=1.2)
    points(NumGenera ~ GlobalBinCenter, globalds, pch=16, col="black", cex=0.8)
    arrows(x0=globalds$GlobalBinCenter, y0=globalds$NumGenera + globalds$NumGRangeThru1, x1=globalds$GlobalBinCenter, y1=globalds$NumGenera, length=0.03, angle=90, code=1, lwd=1.2, col="black")
    
    # number of localities as a bar chart
    for (i in seq(along=unique(globalds$GlobalTimeBins))) { # in global time intervals
      dsi <- globalds[which(globalds$GlobalTimeBins==unique(globalds$GlobalTimeBins)[i]),]
      if (dsi$NumLocs!=0) (polygon(x=c(dsi$GlobalBinStart, dsi$GlobalBinStart, dsi$GlobalBinEnd, dsi$GlobalBinEnd), y=c(-myy.r, dsi$NL_log-myy.r, dsi$NL_log-myy.r, -myy.r), col="lightskyblue", border="navy"))
    }
    for (i in seq(along=unicont)) { # in continental bins
      dsi <- genallrk[which(genallrk$ContinentalBins==unicont[i]),]
      if (dsi$NumLocs!=0) (polygon(x=c(dsi$ContinBinStart, dsi$ContinBinStart, dsi$ContinBinEnd, dsi$ContinBinEnd), y=c(-myy.r, dsi$NL_log-myy.r, dsi$NL_log-myy.r, -myy.r), col=rgb(1,1,1, alpha=0.4)))
    }
    abline(h=-myy.r+log(5)*20, lty=3)
    
    # labelling and legend
    if (j==1) {
      axis(2, at=pretty(genallr$NumGenera), las=1, hadj=0.8, cex.axis=0.9) # y axis with labels for the left sie
      axis(4, at=c(-myy.r, log(c(10,50,100,200))*20-myy.r), labels=FALSE) # y axis without labels for the right side
      text(24.2, 140, LETTERS[2:4][k], pos=4, font=2, cex=1.2) # label multipart figure
    } else { 
      axis(2, at=pretty(genallr$NumGenera), labels=FALSE) # y axis without labels for the left side
      axis(4, at=c(-myy.r, (log(c(10,50,100,200))*20)-myy.r), labels=c(0,10,50,100,200), las=1, hadj=0.2, cex.axis=0.9) # y axis with labels on right side
      text(24.2, 140, LETTERS[6:8][k], pos=4, font=2, cex=1.2) # label multipart figure
    }
    if (k==3) { # for the bottom row of figures, make a nice x axis
      # add boxes and labels showing the stratigraphic stages 
      for (i in seq(along=c(23, mybinends))) {
        if (is.element(i, c(2,4,6,8))) (polygon(x=c(c(23, mybinends)[i-1], c(23, mybinends)[i], c(23, mybinends)[i], c(23, mybinends)[i-1]), y=c(rep(-myy.r-20,2),rep(-myy.r, 2)), col=grey(0.92), border=FALSE, xpd=NA))
      }
      for (i in seq(along=c(23, mybinends))) {
        lines(x=c(c(23, mybinends)[i], c(23, mybinends)[i]), y=c(-myy.r-20, miny-myy.r), xpd=NA)
      }
      lines(x=c(23.8,23.8), y=c(-myy.r,-myy.r-20), xpd=NA)
      lines(x=c(1.8,1.8), y=c(-myy.r,-myy.r-20), xpd=NA)
      lines(x=c(23.8,1.8), y=c(-myy.r-20,-myy.r-20), xpd=NA)
      text(x=rev(sort(unique(genall$GlobalBinCenter[!is.na(genall$GlobalBinCenter)]))), y=-myy.r-10, myshortbins, xpd=NA)
      axis(1, pos=-myy.r-20, padj=-0.5, cex.axis=0.9)
    } else (axis(1, at=seq(25,0,by=-5), labels=FALSE))
    box()
  } # loop through the regions
  if (j==1) { # label the y axes for all panels
    mtext("Number of genera", 2, outer=TRUE, line=1.7, cex=0.7)
  } else {
    mtext("Number of locations", 4, outer=TRUE, line=1.5, cex=0.7)
    # label the x axis for all panels
    mtext(side=1, "Time (million years ago)", outer=TRUE, line=2.7, cex=0.7)
  }
} # loop through the continents
dev.off()


############################################################

# Supplementary figure 2: correlation of diversity estimators
# and correlation of gamma with observed + range-through genera
# plots the figure in two parts that can later be combined

# needs Data_Continents, Data_Regions
# and Output_RangeThruDataContinents.txt, Output_RangeThruDataRegions.txt

setwd("D://Data//Gammadiv//FilesForUpload")

# Part 1: the matrix of correlation graphs for different estimators

# load the datasets and combine them to one
genall <- read.delim("Data_Continents.txt", as.is=TRUE)
names(genall)
genall$Dataset <- "continental"
genall2 <- read.delim("Data_Regions.txt", as.is=TRUE)
names(genall2)
genall2 <- genall2[,-12] # exclude continent column to merge the two
genall2$Dataset <- "regional"
names(genall2) <- names(genall)
genall <- rbind(genall, genall2)
genall <- genall[which(genall$NumLocs>5),]
summary(genall$NumLocs)

# which columns are results for which estimator?
myest <- c(16,18,21,12,14)
mynames <- c("Chao", "Jackknife", "Bootstrap", "Chao1", "ACE")
names(genall)[myest]
summary(genall[,myest])
figlabs <- list(LETTERS[c(1,5,8,11)], LETTERS[c(2,6,9)], LETTERS[c(3,7)], LETTERS[4])

pdf("FigS2-part1.pdf", height=7, width=8)
par(mfcol=c(4,4), mar=c(1,1,0.5,0.5), oma=c(4,4,0,0))
for (i in 1:(length(myest)-1)) {
  myplots <- rev((i+1):length(myest))
  for (k in seq(along=myplots)) {
    # plot the dots for regional and continental datasets across all time intervals
    plot(genall[,myest[i]], genall[,myest[myplots[k]]], type="n", xaxt="n", yaxt="n", xlab="", ylab="", las=1, xlim=c(20,270), ylim=c(20,270))
    points(genall[which(genall$Dataset=="continental"),myest[i]], genall[which(genall$Dataset=="continental"),myest[myplots[k]]], pch=16, cex=1.1)
    points(genall[which(genall$Dataset=="regional"),myest[i]], genall[which(genall$Dataset=="regional"),myest[myplots[k]]], pch=16, col="red3", cex=0.8)
    abline(0,1) # the 1:1 line
    # label the plot, also with the correlation coefficients
    text(15,220, paste("r =", round(cor.test(genall[which(genall$Dataset=="continental"),myest[i]], genall[which(genall$Dataset=="continental"),myest[myplots[k]]])$estimate,3), sep=" "), pos=4)
    text(15,195, paste("r =", round(cor.test(genall[which(genall$Dataset=="regional"),myest[i]], genall[which(genall$Dataset=="regional"),myest[myplots[k]]])$estimate,3), sep=" "), pos=4, col="red3")
    text(25,260, figlabs[[i]][k], font=2, cex=1.2) # label mulit-part figure
    # axes and legends
    if (k==length(myplots)) {
      axis(1)
      mtext(mynames[i], side=1, line=2.6, cex=0.8)
    } else (axis(1, labels=FALSE))
    if (i==1) {
      axis(2, las=1)
      mtext(mynames[myplots[k]], side=2, line=2.8, cex=0.8)
    } else (axis(2, labels=FALSE))
  }
  # make sure the matrix works by putting in empty plots
  if (any(i==2, i==3, i==4)) (plot.new()) 
  if (i==3) {
    # fake a plot to add the legend for red and black dots
    plot(genall[,myest[i]], genall[,myest[myplots[k]]], type="n", xaxt="n", yaxt="n", xlab="", ylab="", bty="n", xlim=c(20,270), ylim=c(20,270))
    legend(19, 160, text.col=c("black", "red3"), legend=c("continental estimates", "regional estimates"), cex=1.2, bty="n")
  }
}
dev.off()

# Part 2: correlation of gamma with observed + range-through genera

# continental data
genall <- read.delim("Data_Continents.txt", as.is=TRUE)
head(genall)
unique(genall$Continent)
genrt <- read.delim("Output_RangeThruDataContinents.txt", as.is=TRUE)
head(genrt)
genall <- merge(genall, genrt, by=c(1:3))
head(genall)
table(genall$Continent, genall$GlobalTimeBins)
table(genall$Continent, genall$ContinentalBins)

# regional data
genallr <- read.delim("Data_Regions.txt", as.is=TRUE)
head(genallr)
unique(genallr$Region)
genrt <- read.delim("Output_RangeThruDataRegions.txt", as.is=TRUE)
head(genrt)
genrt <- genrt[,-length(genrt)] # get rid of continent column for merging
genallr <- merge(genallr, genrt, by=c(1:3), all.x=TRUE)
head(genallr)
table(genallr$Region, genallr$GlobalTimeBins)
table(genallr$Region, genallr$ContinentalBins)
genallr$NumGenera[which(genallr$NumLocs==0)] <- 0

# only use time bins with >5 locations
genallc <- genall[which(genall$NumLocs>5),]
genallc$NumGTotal <- genallc$NumGenera + genallc$NumGRangeThru1
names(genallc)
mygenall <- genallc[,c(1:5,12:22,45,48)] # reduce dataset to columns needed
mygenall$Type <- "Continent"
head(mygenall)
genallreg <- genallr[which(genallr$NumLocs>5),]
genallreg$NumGTotal <- genallreg$NumGenera + genallreg$NumGRangeThru1
names(genallreg)
genallreg <- genallreg[,c(1:5,13:23,46,49)] # reduce dataset to columns needed
genallreg$Type <- "Region"
names(genallreg)[1] <- "Continent"
mygenall <- rbind(mygenall, genallreg)
table(mygenall$Continent, mygenall$GlobalTimeBins)
table(mygenall$Continent, mygenall$ContinentalBins)

summary(mygenall$jack1)
summary(mygenall$S.chao1)
summary(mygenall$NumGTotal)
mygenall[which(mygenall$NumGTotal>mygenall$jack1),]
mygenall[which(mygenall$NumGTotal>mygenall$S.chao1),]
summary(mygenall$NumLocs)
mygenall[which(mygenall$NumLocs < 10),1:3]

pdf("FigS2-part2.pdf", height=2.5, width=7)
par(mfcol=c(1,2), oma=c(4,2.8,1,2.8), mar=c(0.5,1,0.5,1), cex=0.7)
# Jackknife
plot(jack1 ~ NumGTotal, mygenall, type="n", ylim=c(10,270), xlim=c(10,270), ylab="", xlab="", xaxt="n", yaxt="n")
axis(1, padj=-0.5, cex.axis=0.9)
axis(2, las=1, hadj=0.8, cex.axis=0.9)
points(jack1 ~ NumGTotal, mygenall[which(mygenall$Type=="Continent"),], pch=16, cex=1)
points(jack1 ~ NumGTotal, mygenall[which(mygenall$Type=="Region"),], pch=16, col="red3", cex=0.7)
abline(0,1)
# legend(5, 275, text.col=c("black", "red3"), legend=c("continental estimates", "regional estimates"), bty="n")
text(175,65,paste("r =", round(cor.test(mygenall$jack1[mygenall$Type=="Continent"], mygenall$NumGTotal[mygenall$Type=="Continent"])$estimate,3), sep=" "), pos=4)
text(175,45,paste("r =", round(cor.test(mygenall$jack1[mygenall$Type=="Region"], mygenall$NumGTotal[mygenall$Type=="Region"])$estimate,3), sep=" "), pos=4, col="red3")
mtext("Number of observed and range-through genera", side=1, at=300, line=2, cex=0.7)
mtext(expression(paste("Mammalian ", gamma, " diversity")), side=2, outer=TRUE, line=1.7, cex=0.7)
text(0, 265, "L", pos=4, font=2, cex=1.2,xpd=NA)
# Chao1
plot(S.chao1 ~ NumGTotal, mygenall, type="n", ylim=c(10,270), xlim=c(10,270), ylab="", xlab="", xaxt="n", yaxt="n")
axis(1, padj=-0.5, cex.axis=0.9)
points(S.chao1 ~ NumGTotal, mygenall[which(mygenall$Type=="Continent"),], pch=16, cex=1)
points(S.chao1 ~ NumGTotal, mygenall[which(mygenall$Type=="Region"),], pch=16, col="red3", cex=0.7)
abline(0,1)
text(175,65,paste("r = ", round(cor.test(mygenall$S.chao1[mygenall$Type=="Continent"], mygenall$NumGTotal[mygenall$Type=="Continent"])$estimate,3), "0", sep=""), pos=4)
text(175,45,paste("r =", round(cor.test(mygenall$S.chao1[mygenall$Type=="Region"], mygenall$NumGTotal[mygenall$Type=="Region"])$estimate,3), sep=" "), pos=4, col="red3")
text(0, 265, "M", pos=4, font=2, cex=1.2)
dev.off()


#############################################

# Supplementary figure 3: comparison of gamma diversity
# from localities and grid cells

# needs Output_GridCellDataEurope and Data_Regions

setwd("D://Data//Gammadiv//FilesForUpload")
genall <- read.delim("Output_GridCellDataEurope.txt", as.is=TRUE)
head(genall)
# fix a typo... argh.
genall$GlobalTimeBins[which(genall$GlobalTimeBins=="Serravalian")] <- "Serravallian"
gen2 <- read.delim("Data_Regions.txt", as.is=TRUE)
head(gen2)
gen2 <- gen2[,c(1:5,19,20)]
names(gen2)[4:7] <- paste("Or", names(gen2)[4:7], sep=".")
names(gen2)
names(genall)
ds <- merge(genall, gen2, by=c(1:3), all.x=TRUE, all.y=FALSE, sort=FALSE)
head(ds)
# data checking
all(ds$NumLocs==ds$Or.NumLocs)
all(ds$NumGenera[!is.na(ds$NumGenera)]==ds$Or.NumGenera[!is.na(ds$Or.NumGenera)])
my.c <- unique(ds$Region)
# exclude estimates based on too few grid cells or localities
ds[which(ds$NumCells<6),14:24] <- NA
ds[which(ds$Or.NumLocs<6),27:28] <- NA
summary(ds$jack1)
summary(ds$Or.jack1)

# a plot function specially made for the dataset - see Fig 2 code for explanations.
temp.traj <- function(ds, mytime, var.time, var.time.start, var.time.centre, var.time.end, var.y, var.stderr, var.sampling, sampling.min=6, mycol="black", mylwd=1) {
  par(lend="square", lwd=mylwd) # setting some plotting parameters
  for (i in seq(along=mytime)) { # loop through the given time intervals
    dsi <- ds[which(ds[,var.time]==mytime[i]),] # the dataset for interval i
    ynow <- dsi[,var.y] # the value of var.y at time i
    # if the sampling is above the sampling.min threshold and
    # there is a var.y value at time i:
    if (dsi[,var.sampling]>=sampling.min & !is.na(ynow)) {
      # plot the line for time i
      lines(x=c(dsi[,var.time.start], dsi[,var.time.end]), y=rep(ynow, 2), col=mycol)
      # plot the error bar for time i
      arrows(x0=dsi[,var.time.centre], y0=ynow+dsi[,var.stderr], x1=dsi[,var.time.centre], y1=ynow-dsi[,var.stderr], length=0, code=3, col=mycol)
      # if time i is not the first interval
      if (i > 1) {
        # get the var.y value of the previous time interval
        ybefore <- ds[which(ds[var.time]==mytime[i-1]),var.y]
        # if that value was not NA and sampling was good enough, plot a 
        # vertical line to link the value.y of the previous time interval
        # with that of time i
        if (ds[which(ds[,var.time]==mytime[i-1]),var.sampling]>=sampling.min & !is.na(ybefore)) (lines(x=rep(dsi[var.time.start], 2), y=c(ybefore, ynow), col=mycol))
      } # if time i is the first interval, no vertical line is needed
    } # if there is no var.y value at time i, no lines are needed for it
  } # close the loop through time intervals
} # close the function

# make a vector for short bin names to label the x axis
summary(as.factor(ds$GlobalTimeBins))
ds$ShortBin <- substr(ds$GlobalTimeBins, 1, 2) # take the first two letters
ds$ShortBin[ds$ShortBin=="Ea"] <- "EP"
ds$ShortBin[ds$GlobalTimeBins=="Late Pliocene"] <- "LP"
summary(as.factor(ds$ShortBin))
myshortbins <- c(unique(ds$ShortBin[!is.na(ds$ShortBin)])) # make the vector and sort it by temporal order, plus get rid of NA and add Pleistocene
# make a similar vector for end times of bins,
# for labelling the stages into boxes along the x axis
mybinends <- c(unique(ds$GlobalBinEnd[!is.na(ds$GlobalBinEnd)]))

# set plotting parameters for temporal diversity plots
summary(c(ds$jack1+ds$jack1.se, ds$jack1-ds$jack1.se, ds$Or.jack1+ds$Or.jack1.se, ds$Or.jack1-ds$Or.jack1.se))
hist(c(ds$jack1.se, ds$Or.jack1.se))
miny <- 0 # y axis start
maxy <- 210 # y axis end 
unicont <- unique(ds$ContinentalBins[!is.na(ds$ContinentalBins)]) # time bins
# stagger the standard errors by shifting the bars a bit to the right
ds$ContinBinCenter2 <- ds$ContinBinCenter - 0.2
ds$GlobalBinCenter2 <- ds$GlobalBinCenter - 0.2


pdf("FigS3.pdf", width=7, height=5)
par(mfrow=c(2,2), oma=c(4,3.5,1,3.5), mar=c(1.5,0.5,0.5,0.5), lend="square")

# first, the scatterplot of numbers of localities against grid cells, and
# of their gamma diversity from Jackknife

# number of locations against number of grid cells
plot(Or.NumLocs ~ NumCells, ds, xlab="", ylab="", las=1, type="n", cex.axis=0.9, xaxt="n")
points(Or.NumLocs ~ NumCells, ds[which(ds$Region=="Western Europe"),], pch=16, cex=1)
points(Or.NumLocs ~ NumCells, ds[which(ds$Region=="Eastern Europe"),], pch=16, col="red3", cex=1)
axis(1, padj=-0.5, cex.axis=0.9)
mtext("Number of grid cells", side=1, line=2, cex=0.8)
mtext("Number of locations", side=2, line=2.8, cex=0.8)
legend(-5,240, text.col=c("black", "red3"), legend=c("Western Europe", "Eastern Europe"), cex=0.8, bty="n")
text(2, 250, "A", font=2, cex=1.2)
# gamma diversity estimates
plot(Or.jack1 ~ jack1, ds, xlim=c(0,190), ylim=c(0,190), xlab="", ylab="", type="n", xaxt="n", yaxt="n")
points(Or.jack1 ~ jack1, ds[which(ds$Region=="Western Europe"),], pch=16, cex=1)
points(Or.jack1 ~ jack1, ds[which(ds$Region=="Eastern Europe"),], pch=16, col="red3", cex=1)
abline(0,1)
axis(1, padj=-0.5, cex.axis=0.9)
axis(4, cex.axis=0.9, las=1)
mtext(expression(paste(gamma, " diversity from grid cells")), side=1, line=2.2, cex=0.8)
mtext(expression(paste(gamma, " diversity from locations")), side=4, line=2.9, cex=0.8)
text(2, 185, "B", font=2, cex=1.2)

# then, the temporal diversity plots for the 2 regions where the two
# types of diversity (from localities and from grid cells) are compared

par(mar=c(0.5,0.5,2.5,0.5))

# loop through the regions
for (j in seq(along=my.c)) { 
  
  # set up plot with grey bars for every second global time bin
  plot(jack1 ~ GlobalBinCenter, ds, xlim=c(23,1.8), ylim=c(miny, maxy), type="n", xlab="", ylab="", xaxt="n", yaxt="n", xaxs="i", yaxs="i", las=1, bty="n")
  greybins <- c("Aquitanian", "Langhian", "Tortonian", "Early Pliocene")
  for (i in seq(along=greybins)) {
    binstart <- unique(ds$GlobalBinStart[which(!is.na(ds$GlobalBinStart) & ds$GlobalTimeBins==greybins[i])])
    binend <- unique(ds$GlobalBinEnd[which(!is.na(ds$GlobalBinEnd) & ds$GlobalTimeBins==greybins[i])])
    polygon(x=c(binstart, binend, binend, binstart), y=c(rep(0,2),rep(maxy, 2)), col=grey(0.92), border=FALSE)
  }
  
  # extract the data for region j
  genallj <- ds[which(ds$Region==my.c[j]),]
  
  # location Jackknife in orange, with error bars moved slightly
  temp.traj(genallj, mytime=unicont, var.time="ContinentalBins", var.time.start="ContinBinStart", var.time.centre="ContinBinCenter2", var.time.end="ContinBinEnd", var.y="Or.jack1", var.stderr="Or.jack1.se", var.sampling="NumLocs", sampling.min=6, mycol="orange", mylwd=1.2)
  # cell Jackknife in red
  temp.traj(genallj, mytime=unicont, var.time="ContinentalBins", var.time.start="ContinBinStart", var.time.centre="ContinBinCenter", var.time.end="ContinBinEnd", var.y="jack1", var.stderr="jack1.se", var.sampling="NumCells", sampling.min=6, mycol="red3", mylwd=1.2)
  
  # diversity in the global time intervals as a smoothed line
  globalds <- genallj[which(is.na(genallj$ContinentalBins)),]
  # location Jackknife in grey, with error bars moved
  lines(spline(globalds$GlobalBinCenter, y=globalds$Or.jack1, method="natural", n=50), col=grey(0.4), lwd=1.2)
  arrows(x0=globalds$GlobalBinCenter2, y0=globalds$Or.jack1 + globalds$Or.jack1.se, x1=globalds$GlobalBinCenter2, y1=globalds$Or.jack1-globalds$Or.jack1.se, length=0, code=3, lwd=1.2, col=grey(0.4))
  # cell Jackknife in black
  lines(spline(globalds$GlobalBinCenter, y=globalds$jack1, method="natural", n=50), col="black", lwd=1.2)
  arrows(x0=globalds$GlobalBinCenter, y0=globalds$jack1 + globalds$jack1.se, x1=globalds$GlobalBinCenter, y1=globalds$jack1-globalds$jack1.se, length=0, code=3, lwd=1.2, col="black")
  
  # labelling and legend
  par(lwd=1)
  if (j==1) { # add full y axis only for the first plot
    # legend(4, 298, legend=c("stages, Chao1", "stages, Jackknife", "ages/zones, Chao1", "ages/zones, Jackknife"), text.col=c("black", grey(0.4), "red3", "orange"), bty="n", adj=1)
    axis(2, las=1, hadj=0.8, cex.axis=0.9) # y axis with labels for the left column of figures
    text(23, 194, "C", pos=4, font=2, cex=1.2) # label multipart figure
  } else {
    axis(2, labels=FALSE, cex.axis=0.9) # y axis without labels
    text(23, 194, "D", pos=4, font=2, cex=1.2) # label multipart figure
  }
  # add boxes and labels showing the stratigraphic stages 
  for (i in seq(along=c(23, mybinends))) {
    if (is.element(i, c(2,4,6,8,10))) (polygon(x=c(c(23, mybinends)[i-1], c(23, mybinends)[i], c(23, mybinends)[i], c(23, mybinends)[i-1]), y=c(rep(-20,2),rep(miny, 2)), col=grey(0.92), border=FALSE, xpd=NA))
  }
  for (i in seq(along=c(23, mybinends))) {
    lines(x=c(c(23, mybinends)[i], c(23, mybinends)[i]), y=c(-20, miny), xpd=NA)
  }
  lines(x=c(23,1.8), y=c(-20,-20), xpd=NA)
  text(x=c(rev(sort(unique(ds$GlobalBinCenter[!is.na(ds$GlobalBinCenter)])))), y=-10, myshortbins, xpd=NA, cex=0.8)
  axis(1, pos=-20, padj=-0.5, cex.axis=0.9)
  box()
  if (j==1) (mtext(expression(paste("Mammalian ", gamma, " diversity")), 2, line=2.6, cex=0.8)) else { # label the y axis
    mtext(side=1, "Time (million years ago)", line=2.4, outer=TRUE, cex=0.8) # label the x axis  
  }
} # loop through the regions
dev.off()


########################################################################


# Supplementary figure S4: results of the sensitivity analysis for performance
# of Jackknife with present-day data (where true richness is known)

# needs Output_PresentDaySensitivityGammaSummaryLocs and ---Pres

setwd("D://Data//Gammadiv//Results")
ds <- read.delim("Output_PresentDaySensitivityGammaSummaryLocs.txt", as.is=TRUE)
head(ds)
ds$Continent[which(is.na(ds$Continent))] <- "NA" # turn North America into proper data, not NA...
table(ds$Continent, ds$n)
table(ds$Continent, ds$NumLocs)
dsp <- read.delim("Output_PresentDaySensitivityGammaSummaryPres.txt", as.is=TRUE)
head(dsp)
dsp$Continent[which(is.na(dsp$Continent))] <- "NA"
table(dsp$Continent, dsp$n)
table(dsp$Continent, dsp$Preservation)

# get the proportional performance: Jackknife estimate / True genus richness
# for the locations dataset
ds$LocsCont <- paste(ds$Continent, ds$NumLocs, sep="_")
ds$TrueGenera <- 27
ds$TrueGenera[ds$Continent=="NA"] <- 28
ds$PropJack1 <- ds$jack1 / ds$TrueGenera
hist(ds$NumGenera)
hist(ds$jack1)
hist(ds$PropJack1)
# for the preservation dataset
dsp$TrueGenera <- 27
dsp$TrueGenera[dsp$Continent=="NA"] <- 28
dsp$PropJack1 <- dsp$jack1 / dsp$TrueGenera
hist(dsp$PropJack1)
dsp$LocsCont <- paste(dsp$Continent, dsp$Preservation, sep="_")

# plot the median and interquartile range of 100 iterations
# for each NumLocs-by-continent combination:
# make a new data frame just for plotting the locations data
summaries <- unlist(tapply(ds$jack1, ds$LocsCont, quantile))
sumpropies <- unlist(tapply(ds$PropJack1, ds$LocsCont, quantile))
sumnumies <- unlist(tapply(ds$NumGenera, ds$LocsCont, quantile))
nicenames <- c("Jack1", "Prop", "Obs")
for (i in 1:3) { # make a nice new dataframe for the medians, min/max values etc
  if (i==1) (summaries <- summaries) else {
    if (i==2) (summaries <- sumpropies) else (summaries <- sumnumies)
  }
  newdsi <- data.frame(LocsCont=sort(unique(ds$LocsCont)), Min=summaries[seq(1,40,by=5)], Q25=summaries[seq(2,40,by=5)], Median=summaries[seq(3,40,by=5)], Q75=summaries[seq(4,40,by=5)], Max=summaries[seq(5,40,by=5)])
  names(newdsi)[2:6] <- paste(nicenames[i], names(newdsi)[2:6], sep="")
  if (i==1) (newds <- newdsi) else (newds <- merge(newds, newdsi))
}
newds$Continent <- c(rep("Europe",4), rep("North America", 4))
newds$NumLocs <- ds$NumLocs[match(newds$LocsCont, ds$LocsCont)]
newds <- newds[order(newds$Continent, newds$NumLocs),]
newds

# make a new data frame just for plotting the preservation data
summaries <- unlist(tapply(dsp$jack1, dsp$LocsCont, quantile))
sumpropies <- unlist(tapply(dsp$PropJack1, dsp$LocsCont, quantile))
sumnumies <- unlist(tapply(dsp$GeneraPerLocMean, dsp$LocsCont, quantile))
nicenames <- c("Jack1", "Prop", "GenPerLoc")
for (i in 1:3) { # make a nice new dataframe for the medians, min/max values etc
  if (i==1) (summaries <- summaries) else {
    if (i==2) (summaries <- sumpropies) else (summaries <- sumnumies)
  }
  newdsi <- data.frame(LocsCont=sort(unique(dsp$LocsCont)), Min=summaries[seq(1,40,by=5)], Q25=summaries[seq(2,40,by=5)], Median=summaries[seq(3,40,by=5)], Q75=summaries[seq(4,40,by=5)], Max=summaries[seq(5,40,by=5)])
  names(newdsi)[2:6] <- paste(nicenames[i], names(newdsi)[2:6], sep="")
  if (i==1) (newdsp <- newdsi) else (newdsp <- merge(newdsp, newdsi))
}
newdsp$Continent <- c(rep("Europe",4), rep("North America", 4))
newdsp$Preservation <- dsp$Preservation[match(newdsp$LocsCont, dsp$LocsCont)]
newdsp <- newdsp[order(newdsp$Continent, newdsp$Preservation),]
newdsp

pdf("FigS4.pdf", width=7, height=6)
par(mfrow=c(2,2), oma=c(0,2,0,0), mar=c(4.5,2.5,0.5,0.5))
# plot proportional performance (number of estimated genera / true number)
# against number of locations
plot(c(newds$PropMin, newds$PropMax+0.05) ~ c(1:8,1:8), type="n", xlim=c(0.75,8.25), las=1, xlab="", ylab="", xaxt="n")
axis(1, at=c(1:8), labels=newds$NumLocs)
abline(h=1, lwd=1, lty=2, col="black")
arrows(x0=c(5:8,1:4)-0.11, y0=newds$PropMedian, y1=newds$PropMedian, x1=c(5:8,1:4)+0.11, angle=90, length=0, code=3, col=c(rep("red3", 4), rep("black", 4)), lwd=3, lend=2)
arrows(x0=c(5:8,1:4), y0=newds$PropQ25, y1=newds$PropQ75, angle=90, length=0, code=3, col=c(rep("red3", 4), rep("black", 4)), lwd=6, lend=2)
arrows(x0=c(5:8,1:4), y0=newds$PropQ25, y1=newds$PropMin, angle=90, length=0, code=3, col=c(rep("red3", 4), rep("black", 4)))
arrows(x0=c(5:8,1:4), y0=newds$PropQ75, y1=newds$PropMax, angle=90, length=0, code=3, col=c(rep("red3", 4), rep("black", 4)))
text(2.5, 1.17, "North America")
text(6.5, 1.165, "Europe", col="red3")
abline(v=4.5)
mtext("A", side=2, line=3, at=1.21, font=2, las=1)
mtext("Proportion estimated", side=2, line=3, cex=0.9)
mtext("Number of locations", side=1, line=2.5, cex=0.9)
# plot proportional performance (number of estimated genera / true number)
# against preservation proportion for the 100 locations dataset
plot(c(newds$PropMin, newds$PropMax+0.05) ~ c(1:8,1:8), type="n", xlim=c(0.75,8.25), las=1, xlab="", ylab="", xaxt="n", yaxt="n")
axis(1, at=c(1:8), labels=newdsp$Preservation)
axis(2, labels=FALSE)
abline(h=1, lwd=1, lty=2, col="black")
arrows(x0=c(5:8,1:4)-0.11, y0=newdsp$PropMedian, y1=newdsp$PropMedian, x1=c(5:8,1:4)+0.11, angle=90, length=0, code=3, col=c(rep("red3", 4), rep("black", 4)), lwd=3, lend=2)
arrows(x0=c(5:8,1:4), y0=newdsp$PropQ25, y1=newdsp$PropQ75, angle=90, length=0, code=3, col=c(rep("red3", 4), rep("black", 4)), lwd=6, lend=2)
arrows(x0=c(5:8,1:4), y0=newdsp$PropQ25, y1=newdsp$PropMin, angle=90, length=0, code=3, col=c(rep("red3", 4), rep("black", 4)))
arrows(x0=c(5:8,1:4), y0=newdsp$PropQ75, y1=newdsp$PropMax, angle=90, length=0, code=3, col=c(rep("red3", 4), rep("black", 4)))
text(2.5, 1.17, "North America")
text(6.5, 1.165, "Europe", col="red3")
abline(v=4.5)
mtext("B", side=2, line=1, at=1.21, font=2, las=1)
mtext("Preservation (% of genera)", side=1, line=2.5, cex=0.9)
# plot the observed numbers
plot(c(newds$ObsMin, newds$Jack1Max+2) ~ c(1:8,1:8), type="n", las=1, xlim=c(0.75,8.25), xlab="", ylab="", xaxt="n")
axis(1, at=c(1:8), labels=newds$NumLocs)
arrows(x0=c(5:8,1:4)-0.11, y0=newds$ObsMedian, y1=newds$ObsMedian, x1=c(5:8,1:4)+0.11, angle=90, length=0, code=3, col=c(rep("red3", 4), rep("black", 4)), lwd=3, lend=2)
arrows(x0=c(5:8,1:4), y0=newds$ObsQ25, y1=newds$ObsQ75, angle=90, length=0, code=3, col=c(rep("red3", 4), rep("black", 4)), lwd=6, lend=2)
arrows(x0=c(5:8,1:4), y0=newds$ObsQ25, y1=newds$ObsMin, angle=90, length=0, code=3, col=c(rep("red3", 4), rep("black", 4)))
arrows(x0=c(5:8,1:4), y0=newds$ObsQ75, y1=newds$ObsMax, angle=90, length=0, code=3, col=c(rep("red3", 4), rep("black", 4)))
lines(x=c(0.7,4.3), y=c(26,26), col="black", lwd=1, lty=2)
lines(x=c(4.7,8.3), y=c(24,24), col="red3", lwd=1, lty=2)
text(2.5, 33.5, "North America")
text(6.5, 33.3, "Europe", col="red3")
abline(v=4.5)
mtext("C", side=2, line=3, at=34.5, font=2, las=1)
mtext("Number of genera", side=2, line=3, cex=0.9)
mtext("Number of locations", side=1, line=2.5, cex=0.9)
# plot the actual estimates against the number of locations
plot(c(newds$ObsMin, newds$Jack1Max+2) ~ c(1:8,1:8), type="n", xlim=c(0.75,8.25), las=1, xlab="", ylab="", xaxt="n", yaxt="n")
axis(1, at=c(1:8), labels=newds$NumLocs)
axis(2, labels=FALSE)
arrows(x0=c(5:8,1:4)-0.11, y0=newds$Jack1Median, y1=newds$Jack1Median, x1=c(5:8,1:4)+0.11, angle=90, length=0, code=3, col=c(rep("red3", 4), rep("black", 4)), lwd=3, lend=2)
arrows(x0=c(5:8,1:4), y0=newds$Jack1Q25, y1=newds$Jack1Q75, angle=90, length=0, code=3, col=c(rep("red3", 4), rep("black", 4)), lwd=6, lend=2)
arrows(x0=c(5:8,1:4), y0=newds$Jack1Q25, y1=newds$Jack1Min, angle=90, length=0, code=3, col=c(rep("red3", 4), rep("black", 4)))
arrows(x0=c(5:8,1:4), y0=newds$Jack1Q75, y1=newds$Jack1Max, angle=90, length=0, code=3, col=c(rep("red3", 4), rep("black", 4)))
lines(x=c(0.7,4.3), y=c(28,28), col="black", lwd=1, lty=2)
lines(x=c(4.7,8.3), y=c(27,27), col="red3", lwd=1, lty=2)
text(2.5, 33.5, "North America")
text(6.5, 33.3, "Europe", col="red3")
abline(v=4.5)
mtext("D", side=2, line=1, at=34.5, font=2, las=1)
mtext("Number of locations", side=1, line=2.5, cex=0.9)
dev.off()


