
# SUPPLEMENTARY ANALYSES

# line  20 - calculate richness of range-through genera, i.e. those that jump
#           a time bin or more; get some numbers to sum up the results (line 250)
# line 315 - diversity patterns with different location definitions
#           (results line 570)
# line 605 - assessment of Jackknife estimation performance with present day data


# Script by Susanne Fritz, last changed 14th June 2016

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

# exploring the effects of including range-through lineages

# identify lineages that jump a time bin or more
# calculate their richness values for each time bin in continental and
# regional datasets

# needs Data_GenusOccurrence.txt, for exploring results also Data_Continents.txt and Data_Regions.txt
# produces Output_RangeThruDataContinents.txt and Output_RangeThruDataRegions.txt


setwd("D://Data//Gammadiv//FilesForUpload")
ds <- read.delim("Data_GenusOccurrence.txt", as.is=TRUE)
head(ds)
summary(as.factor(ds$Region))
summary(as.factor(ds$GlobalTimeBins))
ds$GlobalTimeBins[which(ds$GlobalTimeBins=="Serravalian")] <- "Serravallian"

# function to make a vector that gives the gaps in the record for a given genus,
# as well as the length of those gaps:
# - vector as long as the vector of abundances for a genus (x) 
# - 0 means the genus occurs in that time interval, or does not occur
# - non-zero entries mark the gaps, and the number gives the length of the gap
get.gaps <- function(x) {
  occ <- as.numeric(x>0) # occurrence as 0-1 vector
  if (length(which(occ>0))>1) { # if genus occurs at all
    lags <- diff(occ) # the difference between vector items in occ
    if (any(lags==-1)) { # if there is any disappearance in the record...
      # ...check that there is also a reappearance
      if (any(lags[which(lags==-1)[1]:length(lags)]==1)) {
        minx <- which(lags==-1)[1] # the first disappearance
        ret <- rep(0, length(1:minx-1)) # return vector up until that first disappearance
        while (length(minx)>0) { # loop until there is no disappearance
          # make sure occurrences between gaps are filled in in the return vector
          if (length(ret)<minx[1]) (ret <- append(ret, rep(0, minx[1] - length(ret))))
          stopi <- which(lags==1)[which(lags==1)>minx[1]] # the end of the gap that starts at minx
          # if there is an end to the gap, append the gap indices to the return vector, else just append 0 until the end of x
          if (length(stopi)>0) (ret <- append(ret, rep(stopi[1] - minx[1], stopi[1]-minx[1]))) else (ret <- append(ret, rep(0, length(x)-length(ret))))
          # get the next disappearance
          if (length(stopi)>0) (minx <- which(lags==-1)[which(lags==-1)>stopi[1]]) else (minx <- numeric(0))
        }
        # make sure ret has right length (this needs to be done if the last gap stops at the last time bin)
        if (length(ret)<length(x)) (ret <- append(ret, rep(0, length(x)-length(ret))))
      } else (ret <- rep(0, length(x)))
    } else (ret <- rep(0, length(x)))
  } else (ret <- rep(0, length(x)))
  return(ret)
}
# testing the function...
get.gaps(c(0,0,0))
get.gaps(c(1,1,1))
get.gaps(c(0,1,1,1,1))
get.gaps(c(0,1,1,1,0))
get.gaps(c(0,0,0,0,1))
get.gaps(c(1,1,1,1,0))
get.gaps(c(0,1,0,0,1))
get.gaps(c(1,1,0,0,1))
get.gaps(c(1,1,1,0,0))
get.gaps(c(1,1,0,1,0))
get.gaps(c(1,0,0,0,1))
get.gaps(c(1,0,1,0,0))
get.gaps(c(1,0,1,0,1))
get.gaps(c(1,0,1,0,0,1))
get.gaps(c(1,0,0,0,1,0))
get.gaps(c(1,0,1,0,1,0,1))
get.gaps(c(1,0,1,0,1,0,1,0))
get.gaps(c(0,1,0,1,0,1,0,1))

# function to count the species with gap occurrences across a genus matrix
# that was made with an apply() of get.gaps()
count.gapsp <- function(x, gaps=1) {
  length(which(x==gaps))
}


# Regions, continental time bins

# get matrix of gap occurrences
unireg <- unique(ds$Region)
ds$ContinentalBins <- as.factor(ds$ContinentalBins)
unicotime <- unique(as.character(ds$ContinentalBins[!is.na(ds$ContinentalBins)]))
unicotime <- unicotime[c(1:18,25,26,19:23,27:30,24,31)]
ds$Genus <- as.factor(ds$Genus)
for (i in seq(along=unireg)) { # loop through the regions
  print(i)
  dsi <- ds[which(ds$Region==unireg[i]),] # data set in region i
  mytabi <- as.matrix(table(dsi$Genus, dsi$ContinentalBins)) # genus-by-time bin table
  mytabi <- mytabi[,match(unicotime, colnames(mytabi))] # resort to right order!
  gapsmati <- apply(mytabi, 1, get.gaps) # get the gaps and how long they are
  rownames(gapsmati) <- paste(unireg[i], colnames(mytabi), sep="_") # keep track of time bins
  if (i==1) (gapsmat <- gapsmati) else (gapsmat <- rbind(gapsmat, gapsmati))
}
str(gapsmat)
output <- expand.grid(unireg, unicotime)
names(output) <- c("Region", "ContinentalBins")
output$Combi <- paste(output$Region, output$ContinentalBins, sep="_")
output[1:21,]
# count the number of species with gaps and fill in the output file
NumRangeThruSp1 <- apply(gapsmat, 1, count.gapsp, gaps=1)
NumRangeThruSp2 <- apply(gapsmat, 1, count.gapsp, gaps=2)
NumRangeThruSp3 <- apply(gapsmat, 1, count.gapsp, gaps=3)
output$Var1 <- NumRangeThruSp1[match(output$Combi, names(NumRangeThruSp1))]
output$Var2 <- NumRangeThruSp2[match(output$Combi, names(NumRangeThruSp2))]
output$Var3 <- NumRangeThruSp3[match(output$Combi, names(NumRangeThruSp3))]
names(output)[which(is.element(names(output), c("Var1", "Var2", "Var3")))] <- paste("NumGRangeThru", c(1:3), sep="")
head(output)
# sort out the continental bins for the right continents
output <- output[,c(1,2,4:6,3)]
output$Combi <- NA # get rid of combi column and instead pretend it's the column for the global time bins
names(output)[6] <- "GlobalTimeBins"
output$Continent <- ds$Continent[match(output$Region, ds$Region)] # add a continent column
# get rid of the wrong bins for each continent (eg the american bins for the european regions)
output$AmBins <- is.element(output$ContinentalBins, as.character(unicotime[17:31]))
summary(as.factor(output$Continent))
summary(output$AmBins)
out1 <- output[which(output$Continent=="North America" & output$AmBins),]
out2 <- output[which(output$Continent!="North America" & !output$AmBins),]
output1 <- rbind(out1, out2)
names(output1)
summary(as.factor(output1$Continent))
output1 <- output1[,c(1,6,2:5,7)]

# Regions, global time bins

# get matrix of gap occurrences
ds$GlobalTimeBins <- as.factor(ds$GlobalTimeBins)
unigltime <- unique(as.character(ds$GlobalTimeBins[!is.na(ds$GlobalTimeBins)]))
for (i in seq(along=unireg)) { # loop through the regions
  print(i)
  dsi <- ds[which(ds$Region==unireg[i]),] # data set in region i
  mytabi <- as.matrix(table(dsi$Genus, dsi$GlobalTimeBins)) # genus-by-time bin table
  mytabi <- mytabi[,match(unigltime, colnames(mytabi))] # resort to right order!
  gapsmati <- apply(mytabi, 1, get.gaps) # get the gaps and how long they are
  rownames(gapsmati) <- paste(unireg[i], colnames(mytabi), sep="_") # keep track of time bins
  if (i==1) (gapsmat <- gapsmati) else (gapsmat <- rbind(gapsmat, gapsmati))
}
str(gapsmat)
output <- expand.grid(unireg, unigltime)
names(output) <- c("Region", "GlobalTimeBins")
output$Combi <- paste(output$Region, output$GlobalTimeBins, sep="_")
output[1:21,]
# count the number of species with gaps and fill in the output file
NumRangeThruSp1 <- apply(gapsmat, 1, count.gapsp, gaps=1)
NumRangeThruSp2 <- apply(gapsmat, 1, count.gapsp, gaps=2)
NumRangeThruSp3 <- apply(gapsmat, 1, count.gapsp, gaps=3)
output$Var1 <- NumRangeThruSp1[match(output$Combi, names(NumRangeThruSp1))]
output$Var2 <- NumRangeThruSp2[match(output$Combi, names(NumRangeThruSp2))]
output$Var3 <- NumRangeThruSp3[match(output$Combi, names(NumRangeThruSp3))]
names(output)[which(is.element(names(output), c("Var1", "Var2", "Var3")))] <- paste("NumGRangeThru", c(1:3), sep="")
head(output)
# sort out the global time dataset to fit the continental time one
output$Combi <- NA # get rid of combi column and instead pretend it's the column for the continental time bins
names(output)[3] <- "ContinentalBins"
output$Continent <- ds$Continent[match(output$Region, ds$Region)] # add a continent column
output <- rbind(output1, output) # combine with the previous output for the continental bins
write.table(output, "Output_RangeThruDataRegions.txt", sep="\t", row.names=FALSE, col.names=TRUE)


# run the same code for the continental data
# Continents, continental time bins

# get matrix of gap occurrences
unireg <- unique(ds$Continent)
ds$ContinentalBins <- as.factor(ds$ContinentalBins)
for (i in seq(along=unireg)) {
  print(i)
  dsi <- ds[which(ds$Continent==unireg[i]),]
  mytabi <- as.matrix(table(dsi$Genus, dsi$ContinentalBins))
  mytabi <- mytabi[,match(unicotime, colnames(mytabi))]
  gapsmati <- apply(mytabi, 1, get.gaps)
  rownames(gapsmati) <- paste(unireg[i], colnames(mytabi), sep="_")
  if (i==1) (gapsmat <- gapsmati) else (gapsmat <- rbind(gapsmat, gapsmati))
}
str(gapsmat)
output <- expand.grid(unireg, unicotime)
names(output) <- c("Continent", "ContinentalBins")
output$Combi <- paste(output$Continent, output$ContinentalBins, sep="_")
output[1:21,]
# count the number of species with gaps and fill in the output file
NumRangeThruSp1 <- apply(gapsmat, 1, count.gapsp, gaps=1)
NumRangeThruSp2 <- apply(gapsmat, 1, count.gapsp, gaps=2)
NumRangeThruSp3 <- apply(gapsmat, 1, count.gapsp, gaps=3)
output$Var1 <- NumRangeThruSp1[match(output$Combi, names(NumRangeThruSp1))]
output$Var2 <- NumRangeThruSp2[match(output$Combi, names(NumRangeThruSp2))]
output$Var3 <- NumRangeThruSp3[match(output$Combi, names(NumRangeThruSp3))]
names(output)[which(is.element(names(output), c("Var1", "Var2", "Var3")))] <- paste("NumGRangeThru", c(1:3), sep="")
head(output)
# sort out the continental bins for the right continents
output <- output[,c(1,3,2,4:6)]
output$Combi <- NA
names(output)[2] <- "GlobalTimeBins"
output$AmBins <- is.element(output$ContinentalBins, as.character(unicotime[17:31]))
summary(as.factor(output$Continent))
summary(output$AmBins)
out1 <- output[which(output$Continent=="North America" & output$AmBins),]
out2 <- output[which(output$Continent!="North America" & !output$AmBins),]
output1 <- rbind(out1, out2)
names(output1)
output1 <- output1[,c(1:6)]

# Continents, global time bins

# get matrix of gap occurrences
for (i in seq(along=unireg)) {
  print(i)
  dsi <- ds[which(ds$Continent==unireg[i]),]
  mytabi <- as.matrix(table(dsi$Genus, dsi$GlobalTimeBins))
  mytabi <- mytabi[,match(unigltime, colnames(mytabi))]
  gapsmati <- apply(mytabi, 1, get.gaps)
  rownames(gapsmati) <- paste(unireg[i], colnames(mytabi), sep="_")
  if (i==1) (gapsmat <- gapsmati) else (gapsmat <- rbind(gapsmat, gapsmati))
}
str(gapsmat)
output <- expand.grid(unireg, unigltime)
names(output) <- c("Continent", "GlobalTimeBins")
output$Combi <- paste(output$Continent, output$GlobalTimeBins, sep="_")
output
# count the number of species with gaps and fill in the output file
NumRangeThruSp1 <- apply(gapsmat, 1, count.gapsp, gaps=1)
NumRangeThruSp2 <- apply(gapsmat, 1, count.gapsp, gaps=2)
NumRangeThruSp3 <- apply(gapsmat, 1, count.gapsp, gaps=3)
output$Var1 <- NumRangeThruSp1[match(output$Combi, names(NumRangeThruSp1))]
output$Var2 <- NumRangeThruSp2[match(output$Combi, names(NumRangeThruSp2))]
output$Var3 <- NumRangeThruSp3[match(output$Combi, names(NumRangeThruSp3))]
names(output)[which(is.element(names(output), c("Var1", "Var2", "Var3")))] <- paste("NumGRangeThru", c(1:3), sep="")
head(output)
# sort out the global time dataset to fit the continental time one
output$Combi <- NA
names(output)[3] <- "ContinentalBins"
output <- rbind(output1, output)
output
write.table(output, "Output_RangeThruDataContinents.txt", sep="\t", row.names=FALSE, col.names=TRUE)


# explore the results

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

# for continents
genallc <- read.delim("Data_Continents.txt", as.is=TRUE)
head(genallc)
genallc[,1:3]
genallc2 <- read.delim("Output_RangeThruDataContinents.txt", as.is=TRUE)
genallc2[,1:3]
table(genallc2$Continent, genallc2$GlobalTimeBins)
table(genallc2$Continent, genallc2$ContinentalBins)
# get number of range-through genera for global and continental intervals
# for single intervals gaps and 2 intervals gaps
summary(genallc2$NumGRangeThru1[!is.na(genallc2$GlobalTimeBins)])
summary(genallc2$NumGRangeThru1[is.na(genallc2$GlobalTimeBins)])
summary(genallc2$NumGRangeThru2[!is.na(genallc2$GlobalTimeBins)])
summary(genallc2$NumGRangeThru2[is.na(genallc2$GlobalTimeBins)])
# explore the differences between observed+range-through and estimated gamma
genallc <- merge(genallc, genallc2, sort=FALSE)
genallc$NumGTotal <- genallc$NumGenera + genallc$NumGRangeThru1
genallc <- genallc[which(genallc$NumLocs>5),]
head(genallc)
summary(genallc$jack1)
summary(genallc$NumGTotal)
cor.test(genallc$jack1, genallc$NumGTotal)
cor.test(genallc$S.chao1, genallc$NumGTotal)

# for regions
genallr <- read.delim("Data_Regions.txt", as.is=TRUE)
head(genallr)
genallc2 <- read.delim("Output_RangeThruDataRegions.txt", as.is=TRUE)
genallc2[,1:3]
table(genallc2$Region, genallc2$GlobalTimeBins)
table(genallc2$Region, genallc2$ContinentalBins)
summary(genallc2$NumGRangeThru1[!is.na(genallc2$GlobalTimeBins)])
summary(genallc2$NumGRangeThru1[is.na(genallc2$GlobalTimeBins)])
summary(genallc2$NumGRangeThru2[!is.na(genallc2$GlobalTimeBins)])
summary(genallc2$NumGRangeThru2[is.na(genallc2$GlobalTimeBins)])
# explore the differences between observed+range-through and estimated gamma
genallr <- merge(genallr, genallc2, sort=FALSE)
genallr$NumGTotal <- genallr$NumGenera + genallr$NumGRangeThru1
rm(genallc2)
genallr <- genallr[which(genallr$NumLocs>5),]
summary(genallr$jack1)
summary(genallr$NumGTotal)
cor.test(genallr$jack1, genallr$NumGTotal)
cor.test(genallr$S.chao1, genallr$NumGTotal)

# check the plots of gamma diversity against the number of observed+range-through
# genera - later incorporated in Fig S2
plot(jack1 ~ NumGTotal, genallc, pch=16)
points(jack1 ~ NumGTotal, genallr, pch=16, col="grey")
abline(0,1)
plot(S.chao1 ~ NumGTotal, genallc, pch=16)
points(S.chao1 ~ NumGTotal, genallr, pch=16, col="grey")
abline(0,1)

# which estimates are below/ close to the 1:1 line?
genallr$DivJ <- genallr$jack1 / genallr$NumGTotal
hist(genallr$DivJ)
genallr[which(genallr$DivJ < 1.1),]
genallr$DivA <- genallr$S.chao1 / genallr$NumGTotal
hist(genallr$DivA)
genallr[which(genallr$DivA < 1),]

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


# Supplementary analysis to assess the effect of different location definitions
# (Europe: locations are actual sites; North America: locations are grouped sites)

# take the locations of Western and Eastern Europe, lump them into the
# one degree grid cells, and rerun the raw diversity and the sampling scripts

# needs: Data_GenusOccurrence (and Data_RegionsGridded.shp for plots)
# produces: Output_GridCellDataEurope

# load the dataset and extract European localities
setwd("D://Data//Gammadiv//FilesForUpload")
ds <- read.delim("Data_GenusOccurrence.txt", as.is=TRUE)
head(ds)
any(duplicated(ds))
summary(ds)
locs <- ds[which(is.element(ds$Region, c("Western Europe", "Eastern Europe"))),]
locs$latlong <- paste(locs$LAT, locs$LONG, sep="_")
length(unique(locs$latlong)) # 1222 locations
length(unique(locs$WM_ID)) # ...into 257 grid cells

# get the grid file for those cells
library(sp)
library(maptools)
library(maps)
mygr <- readShapePoly("Data_RegionsGridded.shp")
plot(mygr)
summary(mygr@data)
all(is.element(locs$WM_ID, mygr$WM_ID))
subgr <- subset(mygr, is.element(mygr$WM_ID, locs$WM_ID)) # the grid subset
pointcols <- rep("red", length(locs[,1])) # make a colour vector for the two regions
pointcols[which(locs$Region=="Western Europe")] <- "blue"

# plot this sampling as a map: which grid cells are filled with localities?
pdf("Fig_LocationGridMapEurope.pdf", width=6, height=5)
plot(mygr, ylim=c(30,60), xlim=c(-10,40), asp=1, bg="lightblue2", border=NA)
map(fill=TRUE, col="white", lty=0, add=TRUE) # background map
map(add=TRUE, col=grey(0.6), interior=FALSE, lwd=0.7) # continent outlines
plot(subgr, col=grey(0.55), border="black", add=TRUE) # the grid cells with data
points(LAT ~ LONG, locs, col=pointcols, pch=16, cex=0.4)
box() # frame the plot
dev.off()

# make a new grid-cell dataset for just Western and Eastern Europe,
# with the information of number of localities and locations
luni <- function(x) (length(unique(x)))
locnums <- tapply(locs$latlong, locs$WM_ID, luni)
hist(locnums)
length(unique(locs$Genus)) # 411 genera in Europe


# first, start makíng the Data_Regions dataset for just W and E Europe,
# with number of grid cells, locations, and genera

ds <- locs
unireg <- unique(ds$Region)
unigltime <- unique(ds$GlobalTimeBins[!is.na(ds$GlobalTimeBins)])
unicotime <- unique(ds$ContinentalBins[!is.na(ds$ContinentalBins)])
ds$ContinentalBins <- as.factor(ds$ContinentalBins)
ds$GlobalTimeBins <- as.factor(ds$GlobalTimeBins)
for (i in seq(along=unireg)) { # loop through the regions
  # global time bins
  dsi <- ds[which(ds$Region==unireg[i]),] # dataset for region i
  numgen <- tapply(dsi$Genus, dsi$GlobalTimeBins, luni) # total number of genera within each global time bin in region i
  numcells <- tapply(dsi$WM_ID, dsi$GlobalTimeBins, luni) # total number of grid cells within each global time bin in region i
  numloc <- tapply(dsi$latlong, dsi$GlobalTimeBins, luni) # total number of locations within each global time bin in region i
  # combine the data together to a new dataset
  newline <- data.frame(Region=rep(unireg[i], length(unigltime)), GlobalTimeBins=unigltime, ContinentalBins=rep(999, length(unigltime)), NumLocs=numloc[match(unigltime, names(numloc))], NumCells=numcells[match(unigltime, names(numcells))], NumGenera=numgen[match(unigltime, names(numgen))])
  # continental time bins
  dsi <- ds[which(ds$Region==unireg[i] & !is.na(ds$ContinentalBins)),] # dataset for region i
  numgen <- tapply(dsi$Genus, dsi$ContinentalBins, luni) # total number of genera within each continental time bin in region i
  numcells <- tapply(dsi$WM_ID, dsi$ContinentalBins, luni) # total number of grid cells within each continental time bin in region i
  numloc <- tapply(dsi$latlong, dsi$ContinentalBins, luni) # total number of locations within each continental time bin in region i
  # combine the data together to a new dataset
  newline2 <- data.frame(Region=rep(unireg[i], length(unicotime)), GlobalTimeBins=rep(999, length(unicotime)), ContinentalBins=unicotime, NumLocs=numloc[match(unicotime, names(numloc))], NumCells=numcells[match(unicotime, names(numcells))], NumGenera=numgen[match(unicotime, names(numgen))])
  newline <- rbind(newline, newline2)
  if (i==1) (genbl <- newline) else (genbl <- rbind(genbl, newline))
}
# new dataset: number of localities and genera for each region and global time bin
summary(genbl)
unique(genbl$Region)
# sort out zeros and NAs
genbl$ContinentalBins[which(genbl$ContinentalBins==999)] <- NA
for (i in 4:length(genbl)) (genbl[which(is.na(genbl[,i])),i] <- 0)
genbl[which(genbl$NumLocs==0),6] <- NA
# append the missing time and region information
genbl$GlobalBinStart <- ds$GlobalBinStart[match(genbl$GlobalTimeBins, ds$GlobalTimeBins)]
genbl$GlobalBinEnd <- ds$GlobalBinEnd[match(genbl$GlobalTimeBins, ds$GlobalTimeBins)]
genbl$GlobalBinCenter <- genbl$GlobalBinStart - (genbl$GlobalBinStart - genbl$GlobalBinEnd)/2
genbl$ContinBinStart <- ds$ContinBinStart[match(genbl$ContinentalBins, ds$ContinentalBins)]
genbl$ContinBinEnd <- ds$ContinBinEnd[match(genbl$ContinentalBins, ds$ContinentalBins)]
genbl$ContinBinCenter <- genbl$ContinBinStart - (genbl$ContinBinStart - genbl$ContinBinEnd)/2
genbl$Continent <- ds$Continent[match(genbl$Region, ds$Region)]
genbl[1:15,]
write.table(genbl, "Output_GridCellDataEurope.txt", row.names=FALSE, col.names=TRUE, sep="\t")


# second, get the abundance and occurrence matrices to run the sampling algorithms on

ds$LocCoTime <- paste(ds$ContinentalBins, ds$WM_ID, sep="_") # unique cell-by-continental-time combination
length(unique(ds$LocCoTime)) # 644 cell-continentalbin-combinations
ds$LocGlTime <- paste(ds$GlobalTimeBins, ds$WM_ID, sep="_") # unique cell-by-global-time combination
length(unique(ds$LocGlTime)) # 503 cell-globalbin-combinations
library(vegan)

# --- regions, global time bins
# a presence-absence matrix for each genus (columns)
# in each cell for each global time bin
is.present <- function(x) {
  if (length(x)>0) (return(1)) else (return(length(x)))
}
ds$Genus <- as.factor(ds$Genus)
str(ds$Genus) # 411 genera
uniloc <- sort(unique(ds$LocGlTime)) # 503 cells
for (i in seq(along=uniloc)) { # loop through the cells
  print(i)
  dsi <- ds[which(ds$LocGlTime==uniloc[i]),] # occurrences for the cell i
  mymati <- t(as.matrix(tapply(dsi$Genus, dsi$Genus, is.present))) # make the matrix for cell i
  rownames(mymati) <- uniloc[i]
  if (i==1) (mymat <- mymati) else (mymat <- rbind(mymat, mymati))
}
str(mymat)
occmatgl <- mymat
# turn this into an abundance matrix for each region
ds$RegionGlTime <- paste(ds$Region, ds$GlobalTimeBins, sep="_")
unireg <- sort(unique(ds$RegionGlTime)) # 16 region-global time bin-combinations
for (i in seq(along=unireg)) {
  # the presence matrix (genus by cell) for just region i
  mati <- mymat[which(is.element(rownames(mymat), ds$LocGlTime[ds$RegionGlTime==unireg[i]])),]
  # get the abundances for just region i
  if (is.matrix(mati)) { # if there is more than one cell
    abmati <- t(as.matrix(colSums(mati, na.rm=TRUE)))
  } else { # if there is only one cell
    abmati <- t(as.matrix(as.numeric(!is.na(mati))))
    colnames(abmati) <- names(mati)
  }
  rownames(abmati) <- unireg[i]
  if (i==1) (abmat <- abmati) else (abmat <- rbind(abmat, abmati))
}
str(abmat)
abmat[1:10,1:20]
amreggl <- abmat

# --- regions, continental time bins
# a presence-absence matrix for each genus (columns)
# in each grid cell for each continental time bin
str(ds$Genus) # 411 genera
uniloc <- sort(unique(ds$LocCoTime)) # 644 cells
for (i in seq(along=uniloc)) { # loop through the cells
  print(i)
  dsi <- ds[which(ds$LocCoTime==uniloc[i]),] # occurrences for the cell i
  mymati <- t(as.matrix(tapply(dsi$Genus, dsi$Genus, is.present))) # make the matrix for cell i
  rownames(mymati) <- uniloc[i]
  if (i==1) (mymat <- mymati) else (mymat <- rbind(mymat, mymati))
}
str(mymat)
head(mymat)
occmatco <- mymat
# turn this into an abundance matrix for each region
ds$RegionCoTime <- paste(ds$Region, ds$ContinentalBins, sep="_")
unireg <- sort(unique(ds$RegionCoTime)) # 31 region - continental time bin - combinations
for (i in seq(along=unireg)) {
  # the presence matrix (genus by cell) for just region i
  mati <- mymat[which(is.element(rownames(mymat), ds$LocCoTime[ds$RegionCoTime==unireg[i]])),]
  # get the abundances for just region i
  if (is.matrix(mati)) {
    abmati <- t(as.matrix(colSums(mati, na.rm=TRUE)))
  } else {
    abmati <- t(as.matrix(as.numeric(!is.na(mati))))
    colnames(abmati) <- names(mati)
  }
  rownames(abmati) <- unireg[i]
  if (i==1) (abmat <- abmati) else (abmat <- rbind(abmat, abmati))
}
str(abmat)
abmat[1:10,1:20]
amregco <- abmat


# third, get the richness estimator data and save them with the regional dataset

regds <- genbl
regds$RegionGlTime <- paste(regds$Region, regds$GlobalTimeBins, sep="_")
regds$RegionCoTime <- paste(regds$Region, regds$ContinentalBins, sep="_")
regds$RowNumber <- seq(along=regds[,1])
head(regds)
zeroize <- function (x) { 
  x[which(is.na(x))] <- 0
  return(x)
}

# GLOBAL TIME BINS
mymat2 <- apply(occmatgl, 2, zeroize) # turn NA into 0 in the presence-absence matrix
mymat2[1:10,1:10]
regpool <- ds$RegionGlTime[match(rownames(mymat2), ds$LocGlTime)] # classification of regions for the presence-absence matrix localities
head(regpool)
# from the occurrence matrix, estimate species pool size for the regions
poolreggl <- specpool(mymat2, pool=regpool)
# from the abundance matrix, estimate species pool size for the regions
estreggl <- t(estimateR(amreggl))
# merge the two datasets
reggl <- merge(estreggl, poolreggl[,-1], by=0, all=TRUE)
row.names(reggl) <- reggl$Row.names
reggl <- reggl[,-1]
# tidy up the the NAs and NaNs
if (any(is.na(reggl[,length(reggl)]))) (reggl[which(is.na(reggl[,length(reggl)])),length(reggl)] <- 0)
head(reggl)
summary(reggl)
# merge with the regional dataset
regds2 <- merge(regds, reggl, by.x="RegionGlTime", by.y=0, all.x=TRUE, all.y=FALSE, sort=FALSE)
head(regds2)

# CONTINENTAL TIME BINS
mymat2 <- apply(occmatco, 2, zeroize) # turn NA into 0 in the presence-absence matrix
mymat2[1:10,1:10]
regpool <- ds$RegionCoTime[match(rownames(mymat2), ds$LocCoTime)] # classification of regions for the presence-absence matrix localities
head(regpool)
# from the occurrence matrix, estimate species pool size for the regions
# for some reason, this does not work for MN01 (only one cell)?!
poolregco <- specpool(mymat2[-1,], pool=regpool[-1])
# from the abundance matrix, estimate species pool size for the regions
estregco <- t(estimateR(amregco))
# merge the two datasets
regco <- merge(estregco, poolregco[,-1], by=0, all=TRUE)
row.names(regco) <- regco$Row.names
regco <- regco[,-1]
# tidy up the the NAs and NaNs
if (any(is.na(regco[,length(regco)]))) (regco[which(is.na(regco[,length(regco)])),length(regco)] <- 0)
summary(regco)

# merge with the regional dataset
for (i in seq(along=regco)) { # loop through the columns because this is a partial match process
  regds2[which(is.element(regds2$RegionCoTime, rownames(regco))),which(names(regds2)==names(regco)[i])] <- regco[match(regds2$RegionCoTime[which(is.element(regds2$RegionCoTime, rownames(regco)))], rownames(regco)),i]
}

# check this worked...
head(regds2)
names(regds2)
reggl["Western Europe_Aquitanian",]
regds2[which(regds2$RegionCoTime=="Western Europe_MN16"),]
regco["Western Europe_MN16",]
all(regds2$NumLocs==regds2$n)
regds2[which(regds2$NumCells!=regds2$n),]
all(regds2$NumGenera[!is.na(regds2$NumGenera)]==regds2$S.obs[!is.na(regds2$S.obs)])

# sort out the dataset and save
regds2 <- regds2[order(regds2$RowNumber),]
names(regds2)
regds2 <- regds2[,c(2:14,18:28)]
head(regds2)
write.table(regds2, "Output_GridCellDataEurope.txt", sep="\t", row.names=FALSE, col.names=TRUE)


# compare the results with grid cells to those with locations

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)])

# 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

# have a look at the differences between Jackknife estimates from cells and locations
summary(ds$jack1)
summary(ds$Or.jack1)
cor.test(ds$jack1, ds$Or.jack1)

ds$Diff <- ds$Or.jack1 - ds$jack1
summary(abs(ds$Diff))
summary(ds$jack1.se)
summary(ds$Or.jack1.se)


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


# Supplementary analysis to assess performance of Jackknife estimation where
# the true genus richness is known: simulate a fossil record from the present day

# We picked all fossil locations (across all continental time intervals) and
# extracted present-day species lists for these. That dataset is uploaded as
# Data_Suppl_PresentDayOccurrence.txt.


# produces Output_Suppl_PresentDaySensitivityGamma.RData,
# Output_PresentDaySensitivityGammaSummaryLocs.txt, Output_PresentDaySensitivityGammaSummaryPres.txt

# CAREFUL - due to the randomizations, if this code is rerun, results may
# slightly differ. The results reported in our manuscript are based on the 
# uploaded output file Output_Suppl_PresentDaySensitivityGamma.RData. The script
# highlights places where using this file would make a difference.


setwd("D://Data//Gammadiv//FilesForUpload")
ds <- read.delim("Data_Suppl_PresentDayOccurrence.txt", as.is=TRUE)
head(ds)
length(unique(ds$Species))
length(unique(ds$Genus)) # 41 genera sampled...
length(unique(ds$Location)) # in 1523 locations (we excluded locations that were only assigned to a stratigraphic stage but could not be assigned to an MN zone or NALMA)
luni <- function(x) (length(unique(x)))
tapply(ds$Location, ds$Continent, luni)
table(ds$Continent, ds$Genus)
length(unique(ds$Genus[which(ds$Continent=="Europe")])) # 24 genera in Europe
length(unique(ds$Genus[which(ds$Continent=="North America")])) # 26 genera in North America

# now, subsample these locations in different ways
# decided to sample 3 sets of locations (10, 50, 100), each 100 times

# make presence-absence matrices for simulated sets
# function to get presence/absence
is.present <- function(x) {
  if (length(x)>0) (return(1)) else (return(0))
}
# define sample sizes and get location lists
samplesizes <- c(10, 50, 100, 200)
locs <- unique(ds$Location) # 1523 locations in total
contlocs <- ds$Continent[match(locs, ds$Location)]
locseu <- locs[which(contlocs=="Europe")]
locsna <- locs[which(contlocs=="North America")]
ds$Genus <- as.factor(ds$Genus)
# run a loop through the sets of locations and make a lot of p-a matrices
for (i in seq(along=samplesizes)) { # loop through the 3 sets
  for (j in 1:100) { # repeat the whole thing 100 times
    print(paste(i, j, sep=" and "))
    locsieu <- locseu[sample(seq(along=locseu), size=samplesizes[i])]
    locsina <- locsna[sample(seq(along=locsna), size=samplesizes[i])]
    for (k in seq(along=locsieu)) { # need to loop through locations
      dsikeu <- ds[which(is.element(ds$Location, locsieu[k])),] # a dataset for just the kth randomly sampled European location in the jth dataset of size i
      # make a p-a matrix for this location
      mymatikeu <- t(as.matrix(tapply(dsikeu$Genus, dsikeu$Genus, is.present)))
      # same for North America
      dsikna <- ds[which(is.element(ds$Location, locsina[k])),]
      mymatikna <- t(as.matrix(tapply(dsikna$Genus, dsikna$Genus, is.present)))
      # unite the two to one matrix
      mymatik <- rbind(mymatikeu, mymatikna)
      rownames(mymatik) <- paste(c("EU", "NA"), j, k, sep="_")
      if (k==1) (mymati <- mymatik) else (mymati <- rbind(mymati, mymatik))
    }
    if (j==1) (mymat <- mymati) else (mymat <- rbind(mymat, mymati))
  }
  if (i==1) (mymat10 <- mymat) else {
    if (i==2) (mymat50 <- mymat) else {
      if (i==3) (mymat100 <- mymat) else (mymat200 <- mymat)
    }
  }
}
str(mymat10)
str(mymat50)
str(mymat100)
str(mymat200)

# Note that these matrices may be different from the data reported in the paper,
# as re-running randomizations introduces variation. 
# Hence, if you want to reproduce the reported results, use the uploaded file:
load("Output_Suppl_PresentDaySensitivityGamma.RData")

# based on the matrices,
# calculate Jackknife estimates for the simulated data 
# and save the dataset
library(vegan)
zeroize <- function (x) { # turn NA into 0 in a matrix
  x[which(is.na(x))] <- 0
  return(x)
}

# loop through the matrices
for (k in seq(along=samplesizes)) { # loop through the 3 sets
  print(paste(samplesizes[k]))
  
  # turn the NAs in the matrix into zeros
  mymat <- apply(get(paste("mymat", samplesizes[k], sep="")), 2, zeroize)

  # get a vector of which locations are in which continent and in which iteration
  mylocs <- rownames(mymat)
  myconts <- character(0)
  for (i in seq(along=mylocs)) {
    myloci <- unlist(strsplit(mylocs[i], "_", fixed=TRUE))
    myconts <- append(myconts, paste(myloci[1], myloci[2], sep="_"))
  }
  print(paste(length(unique(myconts)), "unique continent-by-iteration combinations"))
  
  # calculate Jackknife and process the data into a nice dataframe
  mypool <- specpool(mymat, pool=myconts)
  names(mypool)[1] <- "NumGenera"
  mypool$Continent <- 999
  mypool$Iteration <- 999
  for (i in seq(along=mypool$NumGenera)) {
    iti <- unlist(strsplit(rownames(mypool)[i], "_", fixed=TRUE))
    mypool$Continent[i] <- iti[1]
    mypool$Iteration[i] <- iti[2]
  }
  print(summary(as.factor(mypool$Continent)))
  mypool$NumLocs <- samplesizes[k]
  if (k==1) (ds <- mypool) else (ds <- rbind(ds, mypool))
}
head(ds)
summary(as.factor(ds$Continent))
table(ds$Continent, ds$NumLocs)
write.table(ds, "Output_PresentDaySensitivityGammaSummaryLocs.txt", sep="\t", row.names=FALSE, col.names=TRUE)


# vary the preserved proportion for the dataset with 100 locations
# if the results of the manuscript should be reproduced, this part can be 
# ignored, up to the "load" line 50 lines below.

mymat <- apply(mymat100, 2, zeroize)
mymat[1:10,1:10]
numgen <- rowSums(mymat) # number of genera at each location
hist(numgen) 
summary(numgen)

# function to randomly take out genus occurrences within a locality, given the 
# number to take out from the existing 1s in a 0-1 vector
occ.reduce <- function(x, outx) {
  myout <- sample(seq(along=x[x==1]), size=outx)
  myoutx <- names(x[x==1])[myout]
  x[which(is.element(names(x), myoutx))] <- 0
  return(x)
}

# loop through the  preservation proportions and then the matrix,
# and generate new matrices withdifferent % of genera preserved
pres <- c(0.8, 0.6, 0.4)
for (i in seq(along=pres)) {
  numgensim <- ceiling(numgen * pres[i])
  numgenout <- numgen - numgensim
  mymati <- mymat
  for (k in seq(along=numgenout)) {
    mymati[k,] <- occ.reduce(mymati[k,], numgenout[k])
  }
  if (i==1) (mymatp80 <- mymati) else {
    if (i==2) (mymatp60 <- mymati) else (mymatp40 <- mymati)
  }
}
mymatp100 <- mymat
hist(rowSums(mymatp100))
hist(rowSums(mymatp80))
hist(rowSums(mymatp60))
hist(rowSums(mymatp40))

# CAREFUL!
# the following line overwrites the saved matrices we used for our manuscript
# with those produced above.
save(mymat10, mymat50, mymat100, mymat200, mymatp40, mymatp60, mymatp80, mymatp100, file="Output_PresentDaySensitivityGamma.RData")

# if you want to reproduce the results reported in the paper, simply load this:
load("Output_Suppl_PresentDaySensitivityGamma.RData")

# based on the preservation matrices,
# calculate Jackknife estimates for the simulated data as above for locations
samplesizes <- c(40, 60, 80, 100)
for (k in seq(along=samplesizes)) { # loop through the 4 sets
  print(paste(samplesizes[k]))
  
  # turn the NAs in the matrix into zeros
  mymat <- get(paste("mymatp", samplesizes[k], sep=""))
  numgen <- rowSums(mymat)
  
  # get a vector of which locations are in which continent and in which iteration
  mylocs <- rownames(mymat)
  myconts <- character(0)
  for (i in seq(along=mylocs)) {
    myloci <- unlist(strsplit(mylocs[i], "_", fixed=TRUE))
    myconts <- append(myconts, paste(myloci[1], myloci[2], sep="_"))
  }
  print(paste(length(unique(myconts)), "unique continent-by-iteration combinations"))
  
  # calculate Jackknife and process the data into a nice dataframe
  mypool <- specpool(mymat, pool=myconts)
  names(mypool)[1] <- "NumGenera"
  hist(mypool$NumGenera)
  mypool$Continent <- 999
  mypool$Iteration <- 999
  for (i in seq(along=mypool$NumGenera)) {
    iti <- unlist(strsplit(rownames(mypool)[i], "_", fixed=TRUE))
    mypool$Continent[i] <- iti[1]
    mypool$Iteration[i] <- iti[2]
  }
  print(summary(as.factor(mypool$Continent)))
  mypool$Preservation <- samplesizes[k]
  mypool$GeneraPerLocMin <- tapply(numgen, myconts, min)
  mypool$GeneraPerLocMean <- tapply(numgen, myconts, mean)
  mypool$GeneraPerLocMax <- tapply(numgen, myconts, max)
  if (k==1) (ds <- mypool) else (ds <- rbind(ds, mypool))
}
head(ds)
summary(as.factor(ds$Continent))
table(ds$Continent, ds$Preservation)
table(ds$Continent, ds$n)
write.table(ds, "Output_PresentDaySensitivityGammaSummaryPres.txt", sep="\t", row.names=FALSE, col.names=TRUE)

