
# ANALYSES OF DIVERSITY-PRODUCTIVITY RELATIONSHIP
# based on the Data_Continents and Data_Regions files

# continental (NA, EU) and regional (6 regions)
# with global time bins (stratigraphic stages)
# including supplementary analyses with different gamma diversity estimates
# and different NPP estimates (temporal bin allocation)

# Line  25: continental datasets, prepare data, check them, transform variables
# Line 160: continental datasets, model exploration with best-dates NPP dataset
# Line 475: continental datasets, final models:
#   --   gamma diversity vs NPP (best dates)
#   --   gamma diversity vs NPP (automatic dates)
#   --   gamma diversity vs NPP (best dates) with area accounting for tectonic extension of North America
# Line 550: regional datasets, prepare data, check them, transform variables
# Line 695: regional datasets, model exploration with best-dates NPP dataset
# Line 1085: regional datasets, final models:
#   --   gamma diversity vs NPP (best dates)
#   --   gamma diversity vs NPP (automatic dates)
#   --   gamma diversity vs NPP (best dates) with area accounting for tectonic extension of North America

# produces Output_ModelDataContinents and Output_ModelDataRegions
#   (on which the uploaded Data_AnalysesContinents and Data_AnalysesRegions were based)

# Script by Susanne Fritz, last changed 18th May 2016

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

# CONTINENTS

# check datasets and variable intercorrelations

setwd("D://Data//Gammadiv//FilesForUpload")
ds <- read.delim("Data_Continents.txt", as.is=TRUE)
names(ds)
ds$GlobalBinLength <- ds$GlobalBinStart - ds$GlobalBinEnd
head(ds)
table(ds$Continent, ds$GlobalTimeBins)
plot(Area ~ as.factor(Continent), ds)
plot(Area2 ~ as.factor(Continent), ds)
plot(GlobalBinLength ~ as.factor(Continent), ds)

# correlations of genus richness and number of localities
cor.test(ds$NumGenera, ds$NumLocs)

# correlation of diversity in the different time interval schemes:
# assign the stratigraphic stage diversity to each of its NALMA/MN zones
# first, load the time zone tables
mytimezoneseu <- read.delim("Data_TimezonesEurope.txt", as.is=TRUE)
mytimezoneseu
mytimezonesam <- read.delim("Data_TimezonesAmerica.txt", as.is=TRUE)
mytimezonesam
# make a column that gives the stage for each NALMA and MN zone
ds$match.stage <- -999
ds$match.stage[which(is.element(ds$ContinentalBins, mytimezoneseu$TimeEurasia))] <- mytimezoneseu$GlobalTimeBins[match(ds$ContinentalBins[which(is.element(ds$ContinentalBins, mytimezoneseu$TimeEurasia))], mytimezoneseu$TimeEurasia)]
ds$match.stage[which(is.element(ds$ContinentalBins, mytimezonesam$TimeAmerica2))] <- mytimezonesam$GlobalTimeBins[match(ds$ContinentalBins[which(is.element(ds$ContinentalBins, mytimezonesam$TimeAmerica2))], mytimezonesam$TimeAmerica2)]
ds$match.stage[which(!is.na(ds$GlobalTimeBins))] <- NA
ds[,c(1:3,46)]
summary(as.factor(ds$match.stage))
# make a column with the matching stage-level Jack1 estimate for each continental time bin
ds$jack1.stage <- -999
for (i in seq(along=ds[,1])) { # loop through the dataset
  if (is.na(ds$GlobalTimeBins[i])) { # only do this for lines with continental-time information
    ds$jack1.stage[i] <- ds$jack1[which(ds$Continent==ds$Continent[i] & ds$NumLocs>5 & ds$GlobalTimeBins==ds$match.stage[i])]
  } else (ds$jack1.stage[i] <- NA)
}
summary(ds$jack1.stage)
summary(ds$jack1[which(!is.na(ds$GlobalTimeBins))])
# plot the correlation
plot(jack1 ~ jack1.stage, ds[which(is.na(ds$GlobalTimeBins)),], xlab="gamma in stages", ylab="gamma in ages/zones")
cor.test(ds$jack1[which(is.na(ds$GlobalTimeBins))], ds$jack1.stage[which(is.na(ds$GlobalTimeBins))])
# test the other way around:
# average gamma diversity across the ages/zones within a stage,
# weighted by the length of the age/zone
ds$jack1.conttime <- -999
ds$ContinBinLength <- ds$ContinBinStart - ds$ContinBinEnd
for (i in seq(along=ds[,1])) { # loop through the dataset
  if (!is.na(ds$GlobalTimeBins[i])) { # only do this for lines with global-time information
    dsi <- ds[which(ds$Continent==ds$Continent[i] & ds$NumLocs>5 & is.element(ds$match.stage, ds$GlobalTimeBins[i])),]
    ds$jack1.conttime[i] <- weighted.mean(dsi$jack1, w=dsi$ContinBinLength, na.rm=TRUE)
  } else (ds$jack1.conttime[i] <- NA)
}
summary(ds$jack1.conttime)
summary(ds$jack1[which(is.na(ds$GlobalTimeBins))])
# plot the correlation
plot(jack1 ~ jack1.conttime, ds[which(!is.na(ds$GlobalTimeBins)),], xlab="mean gamma in ages/zones", ylab="gamma in stages")
cor.test(ds$jack1[which(!is.na(ds$GlobalTimeBins))], ds$jack1.conttime[which(!is.na(ds$GlobalTimeBins))])

# get rid of the continental time bins for analyses
ds <- ds[which(!is.na(ds$GlobalTimeBins)),]
table(ds$Continent, ds$GlobalTimeBins)
# get rid of estimates I don't believe
# climate: have to be at least 2 grid cells with data
# mammals: have to be at least 5 localities
ds[which(ds$NumLocs<6),12:22]
ds[which(ds$NumLocs<6),12:22] <- NA
ds[which(ds$NPP_AutoPlants_NumCells<2),25:34]
ds[which(ds$NPP_AutoPlants_NumCells<2),26:34] <- NA
ds[which(ds$NPP_BestPlants_NumCells<2),35:44]

# correlation of the different diversity estimates with each other
cor(ds[,c(5,12,14,16,18,20,21,48),], use="pairwise")
# correlations of the different NPP measures with each other
cor(ds[,c(26,28,36,38),], use="pairwise")
plot(ds$NPP_AutoPlants_meanhdrwght50, ds$NPP_BestPlants_meanhdrwght50)
abline(0,1)

# reduce the dataset to only the columns needed
names(ds)
dats <- ds[,c(1,2,6:8,45,23,24,4,18,19,48,35,38,37,25,28,27)]
dats

# histograms and transformations of the variables

hist(dats$jack1) # gamma diversity - non-normal
shapiro.test(dats$jack1)
shapiro.test(round(dats$jack1, 0))
shapiro.test(log(dats$jack1))
hist(dats$NPP_AutoPlants_meanhdrwght50) # NPP - marginally non-normal for automatic dating
hist(log(dats$NPP_AutoPlants_meanhdrwght50))
shapiro.test(dats$NPP_AutoPlants_meanhdrwght50)
shapiro.test(log(dats$NPP_AutoPlants_meanhdrwght50))
hist(dats$NPP_BestPlants_meanhdrwght50) 
hist(log(dats$NPP_BestPlants_meanhdrwght50))
shapiro.test(dats$NPP_BestPlants_meanhdrwght50) # normal for best dating
shapiro.test(log(dats$NPP_BestPlants_meanhdrwght50))
hist(dats$Area) # area - normal for continents
hist(log(dats$Area))
shapiro.test(dats$Area)
shapiro.test(log(dats$Area))
hist(log(dats$Area2))
shapiro.test(dats$Area2)
shapiro.test(log(dats$Area2))
hist(dats$GlobalBinLength) # length of time bin - normal
hist(log(dats$GlobalBinLength))
shapiro.test(dats$GlobalBinLength)
shapiro.test(log(dats$GlobalBinLength))

dats$Area_lg <- log(dats$Area)
dats$Area2_lg <- log(dats$Area2)
dats$Time_lg <- log(dats$GlobalBinLength)
dats$NPP_AP_mhw50_lg <- log(dats$NPP_AutoPlants_meanhdrwght50)
dats$NPP_BP_mhw50_lg <- log(dats$NPP_BestPlants_meanhdrwght50)
head(dats)

rm(ds) # to make sure I don't use different datasets...


# correlation of explanatory variables

cor(dats[,-c(1:5,9,11,13,15,16,18)], use="pairwise")
# model 1: best plants NPP
cor(dats[,c(10,23,19,21,5)], use="pairwise")
cor(dats[,c(12,23,19,21,5)], use="pairwise")
# model 2: auto plants NPP
cor(dats[,c(10,22,19,21,5)], use="pairwise")
cor(dats[,c(12,22,19,21,5)], use="pairwise")

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

# models for continents
# first, explore different ways of fitting models with
# model 1  gamma diversity vs NPP from plants (best dates)

library(nlme)
dats$Continent <- as.factor(dats$Continent)
dats$jack1.int <- round(dats$jack1, 0)
dats$jack1.ct.int <- round(dats$jack1.conttime, 0)
savedats <- dats
dats <- dats[which(!is.na(dats$NPP_BP_mhw50_lg)),] # dataset with complete cases
summary(dats)


# -- LM and GLM - hugely overfitted, too many variables!!!

m1 <- lm(jack1 ~ NPP_BP_mhw50_lg + Area + GlobalBinCenter + GlobalBinLength + Continent, data=dats)
summary(m1)
par(mfrow=c(2,2))
plot(m1)
m1.1 <- lm(jack1 ~ NPP_BP_mhw50_lg + Area + GlobalBinCenter + Time_lg + Continent, data=dats)
summary(m1.1)
plot(m1.1)
m1.2 <- lm(jack1 ~ NPP_BP_mhw50_lg + Area_lg + GlobalBinCenter + Time_lg + Continent, data=dats)
summary(m1.2)
plot(m1.2)
AIC(m1.2)
m1.3 <- lm(jack1 ~ NPP_BP_mhw50_lg, data=dats)
summary(m1.3)

m2 <- glm(jack1.int ~ NPP_BP_mhw50_lg + Area + GlobalBinCenter + GlobalBinLength + Continent, data=dats, family=poisson)
summary(m2)
plot(m2)
m2.1 <- glm(jack1.int ~ NPP_BP_mhw50_lg + Area + GlobalBinCenter + Time_lg + Continent, data=dats, family=poisson)
summary(m2.1)
plot(m2.1)
m2.2 <- glm(jack1.int ~ NPP_BP_mhw50_lg + Area_lg + GlobalBinCenter + Time_lg + Continent, data=dats, family=poisson)
summary(m2.2)
plot(m2.2)
m2.3 <- glm(jack1.int ~ NPP_BP_mhw50_lg + Area_lg + GlobalBinCenter + Time_lg + Continent, data=dats, family=quasipoisson)
summary(m2.3)
plot(m2.3)
m2.4 <- glm(jack1.int ~ NPP_BP_mhw50_lg + Area_lg + GlobalBinCenter + Time_lg + Continent, data=dats, family=quasipoisson)
summary(m2.4)
plot(m2.4)
m2.5 <- glm(jack1.int ~ NPP_BP_mhw50_lg, data=dats, family=poisson)
summary(m2.5)
m2.6 <- glm(jack1.int ~ NPP_BP_mhw50_lg, data=dats, family=quasipoisson)
summary(m2.6)


# GLMs: trying backwards model selection

m3.1 <- glm(jack1.int ~ Continent * NPP_BP_mhw50_lg + Area_lg + Time_lg + GlobalBinCenter, dats, family=poisson)
summary(m3.1)
m3.2 <- glm(jack1.int ~ Continent * NPP_BP_mhw50_lg + Area_lg + Time_lg + GlobalBinCenter, dats, family=quasipoisson)
summary(m3.2)
m3.3 <- glm(jack1.int ~ Continent + NPP_BP_mhw50_lg + Area_lg + Time_lg + GlobalBinCenter, dats, family=quasipoisson)
summary(m3.3)
anova(m3.2, m3.3, test="F")
m3.4 <- glm(jack1.int ~ Continent + NPP_BP_mhw50_lg + Area_lg + GlobalBinCenter, dats, family=quasipoisson)
summary(m3.4)
anova(m3.3, m3.4, test="F")
m3.5 <- glm(jack1.int ~  NPP_BP_mhw50_lg + Area_lg + GlobalBinCenter, dats, family=quasipoisson)
summary(m3.5)
anova(m3.4, m3.5, test="F")
m3.6 <- glm(jack1.int ~  NPP_BP_mhw50_lg + Area_lg, dats, family=quasipoisson)
summary(m3.6)
anova(m3.5, m3.6, test="F")
m3.7 <- glm(jack1.int ~  NPP_BP_mhw50_lg, dats, family=quasipoisson)
summary(m3.7)
anova(m3.6, m3.7, test="F")
m3.8 <- glm(jack1.int ~  NPP_BP_mhw50_lg + Area_lg + Continent, dats, family=quasipoisson)
summary(m3.8)
anova(m3.6, m3.8, test="F")


# -- GLS models
# much better because variance heterogeneity can be modelled
# NPP, Area, Time as fixed components; test for continent and bin center as random components


# simplest model: basically linear regression
m4 <- gls(jack1.int ~ NPP_BP_mhw50_lg + Area_lg + Time_lg, data = dats, method="ML")
summary(m4)
# look at residuals to check homogeneity of variance
plot(m4)
par(mfrow=c(1,1))
plot(dats$Area_lg, resid(m4))
plot(dats$Continent, resid(m4))
plot(dats$GlobalBinCenter, resid(m4))
plot(dats$Time_lg, resid(m4))
library(lattice)
xyplot(resid(m4) ~ dats$GlobalBinCenter | dats$Continent)

# model controlling for heterogeneity of variance between continents
m4.1 <- gls(jack1.int ~ NPP_BP_mhw50_lg + Area_lg + Time_lg, data = dats, weights=varIdent(form = ~1|Continent), method="ML")
summary(m4.1)
anova(m4, m4.1) # second model not significantly better but AIC only 1.4 higher, 
# so not significantly worse model either.

# temporal autocorrelation: check the residuals for each time series
# and for the multiple time series combinations
dats$MyRes <- residuals(m4, type="normalized")
dats2 <- dats[order(dats$Continent, dats$GlobalBinCenter),]
par(mfrow=c(1,1))
acf(dats2$MyRes[dats2$Continent=="Europe"])
acf(dats2$MyRes[dats2$Continent=="North America"])
mytimeseries <- cbind(dats2$MyRes[dats2$Continent=="Europe"], c(dats2$MyRes[dats2$Continent=="North America"], NA, NA))
acf(mytimeseries, na.action=na.pass)
dats$MyRes <- residuals(m4.1, type="normalized")
dats2 <- dats[order(dats$Continent, dats$GlobalBinCenter),]
acf(dats2$MyRes[dats2$Continent=="Europe"])
acf(dats2$MyRes[dats2$Continent=="North America"])
mytimeseries <- cbind(dats2$MyRes[dats2$Continent=="Europe"], c(dats2$MyRes[dats2$Continent=="North America"], NA, NA))
acf(mytimeseries, na.action=na.pass)

# no serious temporal autocorrelation! Do the model anyway:
m4.2 <- gls(jack1.int ~ NPP_BP_mhw50_lg + Area_lg + Time_lg, data = dats, correlation = corAR1(form = ~ GlobalBinCenter|Continent), weights=varIdent(form = ~1|Continent), method="ML")
summary(m4.2)
plot(m4.2)
anova(m4, m4.1, m4.2)
anova(m4, m4.2)

# reduce fixed components for the best of these 3 models:
# decided for the one with continent because it was not significantly worse than the one
# without, and the residual plots for continents are aweful.
# kick out time because it is non-significant
m4.3 <- gls(jack1.int ~ NPP_BP_mhw50_lg + Area_lg, data = dats, weights=varIdent(form = ~1|Continent), method="ML")
summary(m4.3)
anova(m4.1, m4.3) # no significant change in model when time is kicked out


# with transformed response variable
dats$jack1_lg <- log(dats$jack1.int)

# simple model
m5 <- gls(jack1_lg ~ NPP_BP_mhw50_lg + Area_lg + Time_lg, data = dats, method="ML")
summary(m5)
plot(m5)
plot(dats$Area_lg, resid(m5))
plot(dats$Continent, resid(m5))
plot(dats$GlobalBinCenter, resid(m5))
plot(dats$Time_lg, resid(m5))
xyplot(resid(m5) ~ dats$GlobalBinCenter | dats$Continent)

# model controlling for heterogeneity of variance between continents
m5.1 <- gls(jack1_lg ~ NPP_BP_mhw50_lg + Area_lg + Time_lg, data = dats, weights=varIdent(form = ~1|Continent), method="ML")
summary(m5.1)
anova(m5, m5.1)

# temporal autocorrelation: check the residuals for each time series
# and for the multiple time series combinations
dats$MyRes <- residuals(m5, type="normalized")
dats2 <- dats[order(dats$Continent, dats$GlobalBinCenter),]
acf(dats2$MyRes[dats2$Continent=="Europe"])
acf(dats2$MyRes[dats2$Continent=="North America"])
mytimeseries <- cbind(dats2$MyRes[dats2$Continent=="Europe"], c(dats2$MyRes[dats2$Continent=="North America"], NA, NA))
acf(mytimeseries, na.action=na.pass)
dats$MyRes <- residuals(m5.1, type="normalized")
dats2 <- dats[order(dats$Continent, dats$GlobalBinCenter),]
acf(dats2$MyRes[dats2$Continent=="Europe"])
acf(dats2$MyRes[dats2$Continent=="North America"])
mytimeseries <- cbind(dats2$MyRes[dats2$Continent=="Europe"], c(dats2$MyRes[dats2$Continent=="North America"], NA, NA))
acf(mytimeseries, na.action=na.pass)

# model with temporal autocorrelation
m5.2 <- gls(jack1_lg ~ NPP_BP_mhw50_lg + Area_lg + Time_lg, data = dats, correlation = corAR1(form = ~ GlobalBinCenter|Continent), weights=varIdent(form = ~1|Continent), method="ML")
summary(m5.2)
anova(m5, m5.1, m5.2)
anova(m5, m5.2)

# best model: reduce fixed component
m5.3 <- gls(jack1_lg ~ NPP_BP_mhw50_lg + Area_lg, data = dats, weights=varIdent(form = ~1|Continent), method="ML")
summary(m5.3)
anova(m5.1, m5.3)  # no significant change when kicking out time
m5.4 <- gls(jack1_lg ~ Area_lg, data = dats, weights=varIdent(form = ~1|Continent), method="ML")
summary(m5.4)
anova(m5.3, m5.4) # but model with NPP is significantly better than model without


# -- GLMMs

library(lme4)
library(MuMIn)

# first models without duration of geological stage, varying random and fixed effects
mm0 <- glmer(jack1.int ~ NPP_BP_mhw50_lg + Area_lg + GlobalBinCenter + (1 | Continent), dats, family=poisson)
summary(mm0)
r.squaredGLMM(mm0)
mm1 <- glmer(jack1.int ~ NPP_BP_mhw50_lg + Area_lg + (1 | GlobalBinCenter) + (1 | Continent), dats, family=poisson)
summary(mm1)
r.squaredGLMM(mm1)
anova(mm0, mm1) # probably not legitimate because models differ in random component?!
mm1.1 <- glmer(jack1.int ~ NPP_BP_mhw50_lg + Area_lg + GlobalBinCenter + (1 | GlobalBinCenter) + (1 | Continent), dats, family=poisson)
summary(mm1.1)
r.squaredGLMM(mm1.1)
anova(mm1, mm1.1)
mm2 <- glmer(jack1.int ~ NPP_BP_mhw50_lg + Area_lg + (GlobalBinCenter | Continent), dats, family=poisson)
summary(mm2)
r.squaredGLMM(mm2)
mm3 <- glmer(jack1.int ~ NPP_BP_mhw50_lg + Area_lg + GlobalBinCenter + (GlobalBinCenter | Continent), dats, family=poisson)
summary(mm3)
anova(mm2, mm3)
mm3.1 <- glmer(jack1.int ~ NPP_BP_mhw50_lg + Area_lg + GlobalBinCenter + (1 + GlobalBinCenter | Continent), dats, family=poisson)
summary(mm3.1) # seems to be same model.
mm4 <- glmer(jack1.int ~ NPP_BP_mhw50_lg + Area_lg + Time_lg + GlobalBinCenter + (1 | Continent), dats, family=poisson)
summary(mm4)
r.squaredGLMM(mm4)
anova(mm4, mm0) #compare to model without Time_lg
mm5 <- glmer(jack1.int ~ NPP_BP_mhw50_lg + Area_lg + Time_lg + (1 | GlobalBinCenter) + (1 | Continent), dats, family=poisson)
summary(mm5)
r.squaredGLMM(mm5)
anova(mm5, mm1)
mm6 <- glmer(jack1.int ~ NPP_BP_mhw50_lg + Area_lg + Time_lg + (GlobalBinCenter | Continent), dats, family=poisson)
summary(mm6)
r.squaredGLMM(mm6)
anova(mm6, mm2) # this seems to be the best model of all.
# try out other random component structure: NPP is significant in all cases,
# but the AIC values are higher than in the mm6.
mm7 <- glmer(jack1.int ~ NPP_BP_mhw50_lg + Time_lg + (Area_lg + GlobalBinCenter | Continent), dats, family=poisson)
summary(mm7)
mm7.1 <- glmer(jack1.int ~ NPP_BP_mhw50_lg + (Area_lg + GlobalBinCenter | Continent), dats, family=poisson)
summary(mm7.1)
mm8 <- glmer(jack1.int ~ NPP_BP_mhw50_lg + Area_lg  + (Time_lg + GlobalBinCenter | Continent), dats, family=poisson)
summary(mm8)
mm9 <- glmer(jack1.int ~ NPP_BP_mhw50_lg + (Area_lg + Time_lg + GlobalBinCenter | Continent), dats, family=poisson)
summary(mm9)

# -- repeat GLMMs with the average Jack1 of continental ages/zones
mm0 <- glmer(jack1.ct.int ~ NPP_BP_mhw50_lg + Area_lg + GlobalBinCenter + (1 | Continent), dats, family=poisson)
summary(mm0)
mm1 <- glmer(jack1.ct.int ~ NPP_BP_mhw50_lg + Area_lg + (1 | GlobalBinCenter) + (1 | Continent), dats, family=poisson)
summary(mm1)
mm1.1 <- glmer(jack1.ct.int ~ NPP_BP_mhw50_lg + Area_lg + GlobalBinCenter + (1 | GlobalBinCenter) + (1 | Continent), dats, family=poisson)
summary(mm1.1)
anova(mm1, mm1.1)
mm2 <- glmer(jack1.ct.int ~ NPP_BP_mhw50_lg + Area_lg + (GlobalBinCenter | Continent), dats, family=poisson)
summary(mm2)
mm3 <- glmer(jack1.ct.int ~ NPP_BP_mhw50_lg + Area_lg + GlobalBinCenter + (GlobalBinCenter | Continent), dats, family=poisson)
summary(mm3)
anova(mm2, mm3)
mm3.1 <- glmer(jack1.ct.int ~ NPP_BP_mhw50_lg + Area_lg + GlobalBinCenter + (1 + GlobalBinCenter | Continent), dats, family=poisson)
summary(mm3.1) # seems to be same model.
mm4 <- glmer(jack1.ct.int ~ NPP_BP_mhw50_lg + Area_lg + Time_lg + GlobalBinCenter + (1 | Continent), dats, family=poisson)
summary(mm4)
anova(mm4, mm0) #compare to model without Time_lg
mm5 <- glmer(jack1.ct.int ~ NPP_BP_mhw50_lg + Area_lg + Time_lg + (1 | GlobalBinCenter) + (1 | Continent), dats, family=poisson)
summary(mm5)
anova(mm5, mm1)
mm6 <- glmer(jack1.ct.int ~ NPP_BP_mhw50_lg + Area_lg + Time_lg + (GlobalBinCenter | Continent), dats, family=poisson)
summary(mm6)
anova(mm6, mm2)

# -- GLMMs further: CHECK THESE MODELS with a different R package

library(MCMCglmm)

# set up the chain to store 1000-2000 iterations with the "thin" argument: for nitt=30000 iterations, the default of thin=10 stores 3000 iterations, every 10th one
# default burnin is set to 3000 iterations
mcmc0 <- MCMCglmm(jack1.int ~ NPP_BP_mhw50_lg + Area_lg + GlobalBinCenter, random=~Continent, data=dats, family='poisson', nitt=30000)
summary(mcmc0)
# model checking follows the course notes given in the MCMCglmm library
autocorr(mcmc0$Sol) # check the lowest non-zero lag values are lower than 0.1
autocorr(mcmc0$VCV) # bit high for residual variance (units)
# pdf("ModelComparisonContinents-plots.pdf", height=10, width=6)
plot(mcmc0$Sol) # summary plots for fixed effects
plot(mcmc0$VCV) # summary plots for random effects (and residual variance, "units")
posterior.mode(mcmc0$Sol) # get posterior distribution modes for each fixed effect

mcmc1 <- MCMCglmm(jack1.int ~ NPP_BP_mhw50_lg + Area_lg,  random=~GlobalBinCenter + Continent, data=dats, family='poisson', nitt=30000)
summary(mcmc1)
autocorr(mcmc1$Sol)
autocorr(mcmc1$VCV)
plot(mcmc1$Sol)
plot(mcmc1$VCV)
mcmc1.1 <- MCMCglmm(jack1.int ~ NPP_BP_mhw50_lg + Area_lg + GlobalBinCenter, random=~GlobalBinCenter + Continent, data=dats, family='poisson', nitt=30000)
summary(mcmc1.1)

# this next model does not work with default settings - need other priors... argh.
# dats$Z.BinCenter <- (dats$GlobalBinCenter - mean(dats$GlobalBinCenter)) / sd(dats$GlobalBinCenter)
# mcmc2 <- MCMCglmm(jack1.int ~ NPP_meanhdrwght50_lg + Area_lg,  random=~us(1+Z.BinCenter):Continent, data=dats, family='poisson', nitt=30000)
myprior<-list(R=list(V=1, nu=0), G=list(G1=list(V=diag(2), nu=2)))
mcmc2 <- MCMCglmm(jack1.int ~ NPP_BP_mhw50_lg + Area_lg,  random=~us(1+GlobalBinCenter):Continent, data=dats, family='poisson', nitt=30000, prior=myprior)
summary(mcmc2)
autocorr(mcmc2$Sol)
autocorr(mcmc2$VCV)
plot(mcmc2$Sol)
plot(mcmc2$VCV)

mcmc4 <- MCMCglmm(jack1.int ~ NPP_BP_mhw50_lg + Area_lg + Time_lg + GlobalBinCenter, random=~Continent, data=dats, family='poisson', nitt=30000)
summary(mcmc4)
autocorr(mcmc4$Sol)
autocorr(mcmc4$VCV)
plot(mcmc4$Sol)
plot(mcmc4$VCV)

mcmc5 <- MCMCglmm(jack1.int ~ NPP_BP_mhw50_lg + Area_lg + Time_lg,  random=~GlobalBinCenter + Continent, data=dats, family='poisson', nitt=30000)
summary(mcmc5)
autocorr(mcmc5$Sol)
autocorr(mcmc5$VCV)
plot(mcmc5$Sol)
plot(mcmc5$VCV)

mcmc6 <- MCMCglmm(jack1.int ~ NPP_BP_mhw50_lg + Area_lg + Time_lg,  random=~us(1+GlobalBinCenter):Continent, data=dats, family='poisson', nitt=30000, prior=myprior)
summary(mcmc6)
autocorr(mcmc6$Sol)
autocorr(mcmc6$VCV)
plot(mcmc6$Sol)
plot(mcmc6$VCV)
# dev.off()

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

# final models for continents:
# run models and save residuals

# --   gamma diversity vs NPP from plants (best dates)
mm.bp <- glmer(jack1.int ~ NPP_BP_mhw50_lg + Area_lg + Time_lg + (GlobalBinCenter | Continent), dats, family=poisson)
summary(mm.bp)
r.squaredGLMM(mm.bp)
dats$Mu.BP <- getME(mm.bp, "mu")
head(dats)
datsbp <- dats[,c(1,2,length(dats))]
datsbp$ContBin <- paste(datsbp$Continent, datsbp$GlobalTimeBins, sep="_")

# --   gamma diversity vs NPP from plants (automatic dates)
dats <- savedats
dats <- dats[which(!is.na(dats$NPP_AP_mhw50_lg)),] # dataset with complete cases
summary(dats)
mm.ap <- glmer(jack1.int ~ NPP_AP_mhw50_lg + Area_lg + Time_lg + (GlobalBinCenter | Continent), dats, family=poisson)
summary(mm.ap)
r.squaredGLMM(mm.ap)
dats$Mu.AP <- getME(mm.ap, "mu")
datsap <- dats[,c(1,2,length(dats))]
datsap$ContBin <- paste(datsap$Continent, datsap$GlobalTimeBins, sep="_")

# --   gamma diversity vs NPP from plants (best dates) with area accounting for tectonic extension of North America
dats <- savedats
dats <- dats[which(!is.na(dats$NPP_BP_mhw50_lg)),] # dataset with complete cases
mm.bp2 <- glmer(jack1.int ~ NPP_BP_mhw50_lg + Area2_lg + Time_lg + (GlobalBinCenter | Continent), dats, family=poisson)
summary(mm.bp2)
r.squaredGLMM(mm.bp2)
dats$Mu.BP2 <- getME(mm.bp2, "mu")
datsbp2 <- dats[,c(1,2,length(dats))]
datsbp2$ContBin <- paste(datsbp2$Continent, datsbp2$GlobalTimeBins, sep="_")

# --   gamma diversity as average from continental ages/zones vs NPP from plants (best dates)
dats <- savedats
dats <- dats[which(!is.na(dats$NPP_BP_mhw50_lg)),]
mm.bp <- glmer(jack1.ct.int ~ NPP_BP_mhw50_lg + Area_lg + Time_lg + (GlobalBinCenter | Continent), dats, family=poisson)
summary(mm.bp)
r.squaredGLMM(mm.bp)
dats$Mu.BP.ct <- getME(mm.bp, "mu")
datsbpct <- dats[,c(1,2,length(dats))]
datsbpct$ContBin <- paste(datsbpct$Continent, datsbpct$GlobalTimeBins, sep="_")

# --   gamma diversity as cont. average vs NPP from plants (automatic dates)
dats <- savedats
dats <- dats[which(!is.na(dats$NPP_AP_mhw50_lg)),] # dataset with complete cases
mm.ap <- glmer(jack1.ct.int ~ NPP_AP_mhw50_lg + Area_lg + Time_lg + (GlobalBinCenter | Continent), dats, family=poisson)
summary(mm.ap)
dats$Mu.AP.ct <- getME(mm.ap, "mu")
datsapct <- dats[,c(1,2,length(dats))]
datsapct$ContBin <- paste(datsapct$Continent, datsapct$GlobalTimeBins, sep="_")

# --   gamma diversity as cont. average vs NPP from plants (best dates) with area accounting for tectonic extension of North America
dats <- savedats
dats <- dats[which(!is.na(dats$NPP_BP_mhw50_lg)),] # dataset with complete cases
mm.bp2 <- glmer(jack1.ct.int ~ NPP_BP_mhw50_lg + Area2_lg + Time_lg + (GlobalBinCenter | Continent), dats, family=poisson)
summary(mm.bp2)
dats$Mu.BP2.ct <- getME(mm.bp2, "mu")
dats$ContBin <- paste(dats$Continent, dats$GlobalTimeBins, sep="_")

# save the mu values
savedats$ContBin <- paste(savedats$Continent, savedats$GlobalTimeBins, sep="_")
savedats$Mu.BP <- datsbp$Mu.BP[match(savedats$ContBin, datsbp$ContBin)]
savedats$Mu.BP2 <- datsbp2$Mu.BP2[match(savedats$ContBin, datsbp2$ContBin)]
savedats$Mu.AP <- datsap$Mu.AP[match(savedats$ContBin, datsap$ContBin)]
savedats$Mu.BP.ct <- datsbpct$Mu.BP.ct[match(savedats$ContBin, datsbpct$ContBin)]
savedats$Mu.BP2.ct <- dats$Mu.BP2.ct[match(savedats$ContBin, dats$ContBin)]
savedats$Mu.AP.ct <- datsapct$Mu.AP.ct[match(savedats$ContBin, datsapct$ContBin)]
savedats <- savedats[which(!is.na(savedats$GlobalTimeBins)),]
savedats[,24:length(savedats)]
names(savedats)
write.table(savedats, "Output_ModelDataContinents.txt", sep="\t", row.names=FALSE, col.names=TRUE)


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

# REGIONS

# check datasets and variable intercorrelations

setwd("D://Data//Gammadiv//FilesForUpload")
ds <- read.delim("Data_Regions.txt", as.is=TRUE)
names(ds)
ds$GlobalBinLength <- ds$GlobalBinStart - ds$GlobalBinEnd
head(ds)
table(ds$Region, ds$GlobalTimeBins)
plot(Area ~ as.factor(Region), ds)
plot(Area2 ~ as.factor(Region), ds)
plot(GlobalBinLength ~ as.factor(Region), ds)

# correlations of genus richness and number of localities
cor.test(ds$NumGenera, ds$NumLocs)

# correlation of diversity in the different time interval schemes:
# assign the stratigraphic stage diversity to each of its NALMA/MN zones
# first, load the time zone tables
mytimezoneseu <- read.delim("Data_TimezonesEurope.txt", as.is=TRUE)
mytimezoneseu
mytimezonesam <- read.delim("Data_TimezonesAmerica.txt", as.is=TRUE)
mytimezonesam
# make a column that gives the stage for each NALMA and MN zone
ds$match.stage <- -999
ds$match.stage[which(is.element(ds$ContinentalBins, mytimezoneseu$TimeEurasia))] <- mytimezoneseu$GlobalTimeBins[match(ds$ContinentalBins[which(is.element(ds$ContinentalBins, mytimezoneseu$TimeEurasia))], mytimezoneseu$TimeEurasia)]
ds$match.stage[which(is.element(ds$ContinentalBins, mytimezonesam$TimeAmerica2))] <- mytimezonesam$GlobalTimeBins[match(ds$ContinentalBins[which(is.element(ds$ContinentalBins, mytimezonesam$TimeAmerica2))], mytimezonesam$TimeAmerica2)]
ds$match.stage[which(!is.na(ds$GlobalTimeBins))] <- NA
ds[,c(1:3,47)]
summary(as.factor(ds$match.stage))
# make a column with the matching stage-level Jack1 estimate for each continental time bin
ds$jack1.stage <- -999
for (i in seq(along=ds[,1])) { # loop through the dataset
  if (is.na(ds$GlobalTimeBins[i])) { # only do this for lines with continental-time information
    # only go on if there are continents-time bin data at all
    if (length(which(ds$Region==ds$Region[i] & ds$NumLocs>5 & ds$GlobalTimeBins==ds$match.stage[i]))>0) {
      ds$jack1.stage[i] <- ds$jack1[which(ds$Region==ds$Region[i] & ds$NumLocs>5 & ds$GlobalTimeBins==ds$match.stage[i])]
    } else (ds$jack1.stage[i] <- NA)
  } else (ds$jack1.stage[i] <- NA)
}
summary(ds$jack1.stage)
summary(ds$jack1[which(!is.na(ds$GlobalTimeBins))])
# plot the correlation
plot(jack1 ~ jack1.stage, ds[which(is.na(ds$GlobalTimeBins)),], xlab="gamma in stages", ylab="gamma in ages/zones")
cor.test(ds$jack1[which(is.na(ds$GlobalTimeBins))], ds$jack1.stage[which(is.na(ds$GlobalTimeBins))])
cor.test(ds$jack1[which(is.na(ds$GlobalTimeBins) & is.element(ds$Region, c("Western North America", "Western Europe", "Eastern Europe")))], ds$jack1.stage[which(is.na(ds$GlobalTimeBins) & is.element(ds$Region, c("Western North America", "Western Europe", "Eastern Europe")))])
# test the other way around:
# average gamma diversity across the ages/zones within a stage,
# weighted by the length of the age/zone
ds$jack1.conttime <- -999
ds$ContinBinLength <- ds$ContinBinStart - ds$ContinBinEnd
for (i in seq(along=ds[,1])) { # loop through the dataset
  if (!is.na(ds$GlobalTimeBins[i])) { # only do this for lines with global-time information
    dsi <- ds[which(ds$Region==ds$Region[i] & ds$NumLocs>5 & is.element(ds$match.stage, ds$GlobalTimeBins[i])),]
    if (length(dsi[,1])>0) { # only go on if there are continents-time bin data at all
      ds$jack1.conttime[i] <- weighted.mean(dsi$jack1, w=dsi$ContinBinLength, na.rm=TRUE)
    } else (ds$jack1.conttime[i] <- NA)
  } else (ds$jack1.conttime[i] <- NA)
}
summary(ds$jack1.conttime)
summary(ds$jack1[which(is.na(ds$GlobalTimeBins))])
# plot the correlation
plot(jack1 ~ jack1.conttime, ds[which(!is.na(ds$GlobalTimeBins)),], xlab="mean gamma in ages/zones", ylab="gamma in stages")
cor.test(ds$jack1[which(!is.na(ds$GlobalTimeBins) & is.element(ds$Region, c("Western North America", "Western Europe", "Eastern Europe")))], ds$jack1.conttime[which(!is.na(ds$GlobalTimeBins) & is.element(ds$Region, c("Western North America", "Western Europe", "Eastern Europe")))])

# get rid of the continental time bins for analyses
ds <- ds[which(!is.na(ds$GlobalTimeBins)),]
table(ds$Region, ds$GlobalTimeBins)
# get rid of estimates I don't believe
# climate: have to be at least 2 grid cells with data
# mammals: have to be at least 5 localities
ds[which(ds$NumLocs<6),13:23]
ds[which(ds$NumLocs<6),1:4]
ds[which(ds$NumLocs<6),13:23] <- NA
ds[which(ds$NPP_AutoPlants_NumCells<2),26:35]
ds[which(ds$NPP_AutoPlants_NumCells<2),27:35] <- NA
ds[which(ds$NPP_BestPlants_NumCells<2),36:45]
ds[which(ds$NPP_BestPlants_NumCells<2),37:45] <- NA
# which regions have data after this cleaning?
ds[which(!is.na(ds$NPP_BestPlants_meanhdrwght50) & !is.na(ds$jack1)),1:4]
# reduce dataset to the 3 best regions!
ds <- ds[which(is.element(ds$Region, c("Western Europe", "Eastern Europe", "Western North America"))),]
table(ds$Region, ds$GlobalTimeBins)

# correlation of the different diversity estimates with each other
cor(ds[,c(5,13,15,17,19,21,22,49),], use="pairwise")
# correlations of the different NPP measures with each other
cor(ds[,c(27,29,37,39),], use="pairwise")
plot(ds$NPP_AutoPlants_meanhdrwght50, ds$NPP_BestPlants_meanhdrwght50)
abline(0,1)

# reduce the dataset to only the columns needed
names(ds)
dats <- ds[,c(1,2,6:8,46,12,24,25,4,19,20,49,36,39,38,26,29,28)]
dats

# histograms and transformations of the variables

hist(dats$jack1) # gamma diversity - normal
shapiro.test(dats$jack1)
shapiro.test(round(dats$jack1, 0))
shapiro.test(log(dats$jack1))
hist(dats$NPP_AutoPlants_meanhdrwght50) # NPP - non-normal
hist(log(dats$NPP_AutoPlants_meanhdrwght50))
shapiro.test(dats$NPP_AutoPlants_meanhdrwght50)
shapiro.test(log(dats$NPP_AutoPlants_meanhdrwght50))
hist(dats$NPP_BestPlants_meanhdrwght50)
hist(log(dats$NPP_BestPlants_meanhdrwght50))
shapiro.test(dats$NPP_BestPlants_meanhdrwght50)
shapiro.test(log(dats$NPP_BestPlants_meanhdrwght50))
hist(dats$Area) # area - normal
hist(log(dats$Area))
shapiro.test(dats$Area)
shapiro.test(log(dats$Area))
hist(log(dats$Area2))
shapiro.test(dats$Area2)
shapiro.test(log(dats$Area2))
hist(dats$GlobalBinLength) # length of time bin - normal
hist(log(dats$GlobalBinLength))
shapiro.test(dats$GlobalBinLength)
shapiro.test(log(dats$GlobalBinLength))

dats$Area_lg <- log(dats$Area)
dats$Area2_lg <- log(dats$Area2)
dats$Time_lg <- log(dats$GlobalBinLength)
dats$NPP_AP_mhw50_lg <- log(dats$NPP_AutoPlants_meanhdrwght50)
dats$NPP_BP_mhw50_lg <- log(dats$NPP_BestPlants_meanhdrwght50)
head(dats)

rm(ds) # to make sure I don't use different datasets...

# correlation of explanatory variables

cor(dats[,-c(1:5,7,9,11,12,13,14,16,17,19)], use="pairwise")
# model 1: best plants NPP
cor(dats[,c(11,24,20,22,5)], use="pairwise")
cor(dats[,c(13,24,20,22,5)], use="pairwise")
# model 2: auto plants NPP
cor(dats[,c(11,23,20,22,5)], use="pairwise")
cor(dats[,c(13,23,20,22,5)], use="pairwise")


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

# models for regions


library(nlme)
dats$Region <- as.factor(dats$Region)
dats$Continent <- as.factor(dats$Continent)
dats$jack1.int <- round(dats$jack1, 0)
dats$jack1.ct.int <- round(dats$jack1.conttime, 0)
savedats <- dats
dats <- dats[which(!is.na(dats$NPP_BP_mhw50_lg) & !is.na(dats$jack1)),] # dataset with complete cases
summary(dats)
# Eastern European Burdigalian has no average cont. gamma diversity
dats[which(!is.na(dats$jack1.int) & is.na(dats$jack1.ct.int)),]


# -- LM and GLM - hugely overfitted, too many variables!!!

m1 <- lm(jack1 ~ NPP_BP_mhw50_lg + Area + GlobalBinCenter + GlobalBinLength + Region, data=dats)
summary(m1)
par(mfrow=c(2,2))
plot(m1)
m1.1 <- lm(jack1 ~ NPP_BP_mhw50_lg + Area + GlobalBinCenter + Time_lg + Region, data=dats)
summary(m1.1)
plot(m1.1)
m1.2 <- lm(jack1 ~ NPP_BP_mhw50_lg + Area_lg + GlobalBinCenter + Time_lg + Region, data=dats)
summary(m1.2)
summary.aov(m1.2)
plot(m1.2)
AIC(m1.2)
m1.3 <- lm(jack1 ~ NPP_BP_mhw50_lg, data=dats)
summary(m1.3)
m1.4 <- lm(jack1 ~ NPP_BP_mhw50_lg * Region, data=dats)
summary.aov(m1.4)
m1.5 <- lm(jack1 ~ NPP_BP_mhw50_lg + Area_lg, data=dats)
summary(m1.5)

m2 <- glm(jack1.int ~ NPP_BP_mhw50_lg + Area + GlobalBinCenter + GlobalBinLength + Region, data=dats, family=poisson)
summary(m2)
plot(m2)
m2.1 <- glm(jack1.int ~ NPP_BP_mhw50_lg + Area + GlobalBinCenter + Time_lg + Region, data=dats, family=poisson)
summary(m2.1)
plot(m2.1)
m2.2 <- glm(jack1.int ~ NPP_BP_mhw50_lg + Area_lg + GlobalBinCenter + Time_lg + Region, data=dats, family=poisson)
summary(m2.2)
plot(m2.2)
m2.3 <- glm(jack1.int ~ NPP_BP_mhw50_lg + Area + GlobalBinCenter + Time_lg + Region, data=dats, family=quasipoisson)
summary(m2.3)
plot(m2.3)
m2.4 <- glm(jack1.int ~ NPP_BP_mhw50_lg + Area_lg + GlobalBinCenter + Time_lg + Region, data=dats, family=quasipoisson)
summary(m2.4)
AIC(m2.4)
plot(m2.4)

# GLMs: trying backwards model selection

m3.1 <- glm(jack1.int ~ Region * NPP_BP_mhw50_lg + Area_lg + Time_lg + GlobalBinCenter, dats, family=poisson)
summary(m3.1)
m3.2 <- glm(jack1.int ~ Region * NPP_BP_mhw50_lg + Area_lg + Time_lg + GlobalBinCenter, dats, family=quasipoisson)
summary(m3.2)
m3.3 <- glm(jack1.int ~ Region + NPP_BP_mhw50_lg + Area_lg + Time_lg + GlobalBinCenter, dats, family=quasipoisson)
summary(m3.3)
anova(m3.2, m3.3, test="F")
m3.4 <- glm(jack1.int ~ Region + NPP_BP_mhw50_lg + Area_lg + GlobalBinCenter, dats, family=quasipoisson)
summary(m3.4)
anova(m3.3, m3.4, test="F")
m3.5 <- glm(jack1.int ~  NPP_BP_mhw50_lg + Area_lg + GlobalBinCenter, dats, family=quasipoisson)
summary(m3.5)
anova(m3.4, m3.5, test="F")
m3.6 <- glm(jack1.int ~  NPP_BP_mhw50_lg + Area_lg, dats, family=quasipoisson)
summary(m3.6)
anova(m3.5, m3.6, test="F")
m3.7 <- glm(jack1.int ~  NPP_BP_mhw50_lg, dats, family=quasipoisson)
summary(m3.7)
anova(m3.6, m3.7, test="F")
m3.8 <- glm(jack1.int ~  Area_lg, dats, family=quasipoisson)
summary(m3.8)
anova(m3.6, m3.8, test="F")
m3.9 <- glm(jack1.int ~  NPP_BP_mhw50_lg + Area_lg + Region, dats, family=quasipoisson)
summary(m3.9)
anova(m3.6, m3.9, test="F")


# -- GLS models
# much better because variance heterogeneity can be controlled for
# NPP, Area, Time as fixed components; test for continent and bin center as random components


# simplest model: basically linear regression
m4 <- gls(jack1.int ~ NPP_BP_mhw50_lg + Area_lg + Time_lg, data = dats, method="ML")
summary(m4)
# look at residuals to check homogeneity of variance
plot(m4)
par(mfrow=c(1,1))
plot(dats$Area_lg, resid(m4))
plot(dats$Region, resid(m4))
plot(dats$GlobalBinCenter, resid(m4))
plot(dats$Time_lg, resid(m4))
library(lattice)
xyplot(resid(m4) ~ dats$GlobalBinCenter | dats$Region)

# model controlling for heterogeneity of variance between continents
m4.1 <- gls(jack1.int ~ NPP_BP_mhw50_lg + Area_lg + Time_lg, data = dats, weights=varIdent(form = ~1|Continent), method="ML")
summary(m4.1)
anova(m4, m4.1) # second model marignally significantly better!
# add the nested region / Continent structure - I can't include grouping factors though!
m4.2 <- gls(jack1.int ~ NPP_BP_mhw50_lg + Area_lg + Time_lg, data = dats, weights=varIdent(form = ~1|Continent|Region), method="ML")
summary(m4.2)
anova(m4, m4.2)
anova(m4.1, m4.2)
# effectively the same - just regions...
m4.3 <- gls(jack1.int ~ NPP_BP_mhw50_lg + Area_lg + Time_lg, data = dats, weights=varIdent(form = ~1|Region), method="ML")
summary(m4.3)
anova(m4, m4.3)
anova(m4.1, m4.3)

# temporal autocorrelation: check the residuals for each time series
# and for the multiple time series combinations
dats$MyRes <- residuals(m4, type="normalized")
dats2 <- dats[order(dats$Region, dats$GlobalBinCenter),]
acf(dats2$MyRes[dats2$Region=="Western Europe"])
acf(dats2$MyRes[dats2$Region=="Eastern Europe"])
acf(dats2$MyRes[dats2$Region=="Western North America"])
mytimeseries <- cbind(dats2$MyRes[dats2$Region=="Western Europe"], c(NA, dats2$MyRes[dats2$Region=="Eastern Europe"], NA), c(dats2$MyRes[dats2$Region=="Western North America"], NA, NA))
acf(mytimeseries, na.action=na.pass)
dats$MyRes <- residuals(m4.3, type="normalized")
dats2 <- dats[order(dats$Region, dats$GlobalBinCenter),]
acf(dats2$MyRes[dats2$Region=="Western Europe"])
acf(dats2$MyRes[dats2$Region=="Eastern Europe"])
acf(dats2$MyRes[dats2$Region=="Western North America"])
mytimeseries <- cbind(dats2$MyRes[dats2$Region=="Western Europe"], c(NA, dats2$MyRes[dats2$Region=="Eastern Europe"], NA), c(dats2$MyRes[dats2$Region=="Western North America"], NA, NA))
acf(mytimeseries, na.action=na.pass)

# no significant temporal autocorrelation! Do the model anyway:
m4.4 <- gls(jack1.int ~ NPP_BP_mhw50_lg + Area_lg + Time_lg, data = dats, correlation = corAR1(form = ~ GlobalBinCenter|Region), weights=varIdent(form = ~1|Region), method="ML")
summary(m4.4)
plot(m4.4)
anova(m4, m4.3, m4.4)
anova(m4, m4.4)

# reduce fixed components for the best of these 3 models:
# kick out time because it is non-significant
m4.5 <- gls(jack1.int ~ NPP_BP_mhw50_lg + Area_lg, data = dats, weights=varIdent(form = ~1|Region), method="ML")
summary(m4.5)
anova(m4.3, m4.5)
# actually also NPP is non-significant...
m4.5a <- gls(jack1.int ~ Area_lg, data = dats, weights=varIdent(form = ~1|Region), method="ML")
summary(m4.5a)
anova(m4.5, m4.5a)
# check the same is true if I put the continent structure in, not the regions
m4.6 <- gls(jack1.int ~ NPP_BP_mhw50_lg + Area_lg, data = dats, weights=varIdent(form = ~1|Continent), method="ML")
summary(m4.6)
anova(m4.1, m4.6) 
# check the best model with temporal autocorrelation
m4.7 <- gls(jack1.int ~ NPP_BP_mhw50_lg + Area_lg, data = dats, correlation = corAR1(form = ~ GlobalBinCenter|Region), weights=varIdent(form = ~1|Region), method="ML")
summary(m4.7)
anova(m4.4, m4.7)
m4.8 <- gls(jack1.int ~ Area_lg, data = dats, correlation = corAR1(form = ~ GlobalBinCenter|Region), weights=varIdent(form = ~1|Region), method="ML")
summary(m4.8)
anova(m4.7, m4.8)

# with transformed response variable
dats$jack1_lg <- log(dats$jack1.int)

# simple model
m5 <- gls(jack1_lg ~ NPP_BP_mhw50_lg + Area_lg + Time_lg, data = dats, method="ML")
summary(m5)
plot(m5)
plot(dats$Area_lg, resid(m5))
plot(dats$Continent, resid(m5))
plot(dats$Region, resid(m5))
plot(dats$GlobalBinCenter, resid(m5))
plot(dats$Time_lg, resid(m5))
xyplot(resid(m5) ~ dats$GlobalBinCenter | dats$Region)

# model controlling for heterogeneity of variance between continents
m5.1 <- gls(jack1_lg ~ NPP_BP_mhw50_lg + Area_lg + Time_lg, data = dats, weights=varIdent(form = ~1|Continent), method="ML")
summary(m5.1)
anova(m5, m5.1)
# model controlling for heterogeneity of variance between regions
m5.2 <- gls(jack1_lg ~ NPP_BP_mhw50_lg + Area_lg + Time_lg, data = dats, weights=varIdent(form = ~1|Region), method="ML")
summary(m5.2)
anova(m5, m5.2) # second model significantly better!

# temporal autocorrelation: check the residuals for each time series
# and for the multiple time series combinations
dats$MyRes <- residuals(m5, type="normalized")
dats2 <- dats[order(dats$Region, dats$GlobalBinCenter),]
acf(dats2$MyRes[dats2$Region=="Western Europe"])
acf(dats2$MyRes[dats2$Region=="Eastern Europe"])
acf(dats2$MyRes[dats2$Region=="Western North America"])
mytimeseries <- cbind(dats2$MyRes[dats2$Region=="Western Europe"], c(NA, dats2$MyRes[dats2$Region=="Eastern Europe"], NA), c(dats2$MyRes[dats2$Region=="Western North America"], NA, NA))
acf(mytimeseries, na.action=na.pass)
dats$MyRes <- residuals(m5.2, type="normalized")
dats2 <- dats[order(dats$Region, dats$GlobalBinCenter),]
acf(dats2$MyRes[dats2$Region=="Western Europe"])
acf(dats2$MyRes[dats2$Region=="Eastern Europe"])
acf(dats2$MyRes[dats2$Region=="Western North America"])
mytimeseries <- cbind(dats2$MyRes[dats2$Region=="Western Europe"], c(NA, dats2$MyRes[dats2$Region=="Eastern Europe"], NA), c(dats2$MyRes[dats2$Region=="Western North America"], NA, NA))
acf(mytimeseries, na.action=na.pass)

# models with temporal autocorrelation: use only the first lag
m5.3 <- gls(jack1_lg ~ NPP_BP_mhw50_lg + Area_lg + Time_lg, data = dats, correlation = corAR1(form = ~ GlobalBinCenter|Region), weights=varIdent(form = ~1|Region), method="ML")
summary(m5.3)
anova(m5, m5.2, m5.3)
anova(m5, m5.3)

# best model: reduce fixed component
m5.4 <- gls(jack1_lg ~ NPP_BP_mhw50_lg + Area_lg, data = dats, weights=varIdent(form = ~1|Region), method="ML")
summary(m5.4)
anova(m5.2, m5.4) 
m5.5 <- gls(jack1_lg ~ Area_lg, data = dats, weights=varIdent(form = ~1|Region), method="ML")
summary(m5.5)
anova(m5.4, m5.5) 


# -- GLMMs

library(lme4)
library(MuMIn)

# first models without duration of geological stage, varying random and fixed effects
mm0 <- glmer(jack1.int ~ NPP_BP_mhw50_lg + Area_lg + GlobalBinCenter + (1 | Region), dats, family=poisson)
summary(mm0)
r.squaredGLMM(mm0)
# try out the nested continent/region structure as the random component
mm0a <- glmer(jack1.int ~ NPP_BP_mhw50_lg + Area_lg + GlobalBinCenter + (1 | Continent/Region), dats, family=poisson)
summary(mm0a)

mm1 <- glmer(jack1.int ~ NPP_BP_mhw50_lg + Area_lg + (1 | GlobalBinCenter) + (1 | Region), dats, family=poisson)
summary(mm1)
r.squaredGLMM(mm1)
mm1a <- glmer(jack1.int ~ NPP_BP_mhw50_lg + Area_lg + (1 | GlobalBinCenter) + (1 | Continent/Region), dats, family=poisson)
summary(mm1a)
mm1.1 <- glmer(jack1.int ~ NPP_BP_mhw50_lg + Area_lg + GlobalBinCenter + (1 | GlobalBinCenter) + (1 | Region), dats, family=poisson)
summary(mm1.1)
anova(mm1, mm1.1)
mm2 <- glmer(jack1.int ~ NPP_BP_mhw50_lg + Area_lg + (GlobalBinCenter | Region), dats, family=poisson)
summary(mm2)
r.squaredGLMM(mm2)
mm2a <- glmer(jack1.int ~ NPP_BP_mhw50_lg + Area_lg + (GlobalBinCenter | Continent/Region), dats, family=poisson)
summary(mm2a)
mm3 <- glmer(jack1.int ~ NPP_BP_mhw50_lg + Area_lg + GlobalBinCenter + (GlobalBinCenter | Region), dats, family=poisson)
summary(mm3)
anova(mm2, mm3)

mm4 <- glmer(jack1.int ~ NPP_BP_mhw50_lg + Area_lg + Time_lg + GlobalBinCenter + (1 | Region), dats, family=poisson)
summary(mm4)
r.squaredGLMM(mm4)
anova(mm4, mm0) #compare to model without Time_lg - significantly better.
mm4a <- glmer(jack1.int ~ NPP_BP_mhw50_lg + Area_lg + Time_lg + GlobalBinCenter + (1 | Continent/Region), dats, family=poisson)
summary(mm4a)
mm5 <- glmer(jack1.int ~ NPP_BP_mhw50_lg + Area_lg + Time_lg + (1 | GlobalBinCenter) + (1 | Region), dats, family=poisson)
summary(mm5)
r.squaredGLMM(mm5)
anova(mm5, mm1)
mm5a <- glmer(jack1.int ~ NPP_BP_mhw50_lg + Area_lg + Time_lg + (1 | GlobalBinCenter) + (1 | Continent/Region), dats, family=poisson)
summary(mm5a)
mm6 <- glmer(jack1.int ~ NPP_BP_mhw50_lg + Area_lg + Time_lg + (GlobalBinCenter | Region), dats, family=poisson)
summary(mm6)
r.squaredGLMM(mm6)
anova(mm6, mm2) # this seems to be the best model of all.
mm6a <- glmer(jack1.int ~ NPP_BP_mhw50_lg + Area_lg + Time_lg + (GlobalBinCenter | Continent/Region), dats, family=poisson)
summary(mm6a)
# try out other random component structure: NPP is significant in mm7 and mm8 only.
mm7 <- glmer(jack1.int ~ NPP_BP_mhw50_lg + Time_lg + (Area_lg + GlobalBinCenter | Region), dats, family=poisson)
summary(mm7)
mm7.1 <- glmer(jack1.int ~ NPP_BP_mhw50_lg + (Area_lg + GlobalBinCenter | Region), dats, family=poisson)
summary(mm7.1)
mm8 <- glmer(jack1.int ~ NPP_BP_mhw50_lg + Area_lg  + (Time_lg + GlobalBinCenter | Region), dats, family=poisson)
summary(mm8)
mm9 <- glmer(jack1.int ~ NPP_BP_mhw50_lg + (Area_lg + Time_lg + GlobalBinCenter | Region), dats, family=poisson)
summary(mm9)

# check the best models

dats$MyRes <- residuals(mm6, type="response")
hist(dats$MyRes)
boxplot(MyRes ~ Region, data=dats, las=2, cex.axis=0.5)
dats$MyFitted <- fitted(mm6)
coplot(MyRes ~ MyFitted | Region, data=dats)
coplot(MyRes ~ NPP_BP_mhw50_lg | Region, data=dats)
coplot(MyRes ~ Area_lg | Region, data=dats)
coplot(MyRes ~ Time_lg | Region, data=dats)
coplot(MyRes ~ GlobalBinCenter | Region, data=dats)

dats$MyRes <- residuals(mm2, type="response")
hist(dats$MyRes)
boxplot(MyRes ~ Region, data=dats, las=2, cex.axis=0.5)
dats$MyFitted <- fitted(mm2)
coplot(MyRes ~ MyFitted | Region, data=dats)
coplot(MyRes ~ NPP_BP_mhw50_lg | Region, data=dats)
coplot(MyRes ~ Area_lg | Region, data=dats)
coplot(MyRes ~ Time_lg | Region, data=dats)
coplot(MyRes ~ GlobalBinCenter | Region, data=dats)


# -- repeat GLMMs with the average Jack1 of continental ages/zones
mm0 <- glmer(jack1.ct.int ~ NPP_BP_mhw50_lg + Area_lg + GlobalBinCenter + (1 | Region), dats, family=poisson)
summary(mm0)
mm1 <- glmer(jack1.ct.int ~ NPP_BP_mhw50_lg + Area_lg + (1 | GlobalBinCenter) + (1 | Region), dats, family=poisson)
summary(mm1)
mm1a <- glmer(jack1.ct.int ~ NPP_BP_mhw50_lg + Area_lg + (1 | GlobalBinCenter) + (1 | Continent/Region), dats, family=poisson)
summary(mm1a)
mm1.1 <- glmer(jack1.ct.int ~ NPP_BP_mhw50_lg + Area_lg + GlobalBinCenter + (1 | GlobalBinCenter) + (1 | Region), dats, family=poisson)
summary(mm1.1)
anova(mm1, mm1.1)
mm2 <- glmer(jack1.ct.int ~ NPP_BP_mhw50_lg + Area_lg + (GlobalBinCenter | Region), dats, family=poisson)
summary(mm2)
mm2a <- glmer(jack1.ct.int ~ NPP_BP_mhw50_lg + Area_lg + (GlobalBinCenter | Continent/Region), dats, family=poisson)
summary(mm2a)
mm3 <- glmer(jack1.ct.int ~ NPP_BP_mhw50_lg + Area_lg + GlobalBinCenter + (GlobalBinCenter | Region), dats, family=poisson)
summary(mm3)
anova(mm2, mm3)
mm4 <- glmer(jack1.ct.int ~ NPP_BP_mhw50_lg + Area_lg + Time_lg + GlobalBinCenter + (1 | Region), dats, family=poisson)
summary(mm4)
anova(mm4, mm0) #compare to model without Time_lg - significantly better.
mm4a <- glmer(jack1.ct.int ~ NPP_BP_mhw50_lg + Area_lg + Time_lg + GlobalBinCenter + (1 | Continent/Region), dats, family=poisson)
summary(mm4a)
mm5 <- glmer(jack1.ct.int ~ NPP_BP_mhw50_lg + Area_lg + Time_lg + (1 | GlobalBinCenter) + (1 | Region), dats, family=poisson)
summary(mm5)
anova(mm5, mm1)
mm5a <- glmer(jack1.ct.int ~ NPP_BP_mhw50_lg + Area_lg + Time_lg + (1 | GlobalBinCenter) + (1 | Continent/Region), dats, family=poisson)
summary(mm5a)
mm6 <- glmer(jack1.ct.int ~ NPP_BP_mhw50_lg + Area_lg + Time_lg + (GlobalBinCenter | Region), dats, family=poisson)
summary(mm6)
anova(mm6, mm2) # this seems to be the best model of all.
mm6a <- glmer(jack1.ct.int ~ NPP_BP_mhw50_lg + Area_lg + Time_lg + (GlobalBinCenter | Continent/Region), dats, family=poisson)
summary(mm6a)


# -- GLMMs further: CHECK THESE MODELS with a different R package

library(MCMCglmm)

# set up the chain to store 1000-2000 iterations with the "thin" argument: for nitt=30000 iterations, the default of thin=10 stores 3000 iterations, every 10th one
# default burnin is set to 3000 iterations
mcmc0 <- MCMCglmm(jack1.int ~ NPP_BP_mhw50_lg + Area_lg + GlobalBinCenter, random=~Region, data=dats, family='poisson', nitt=30000)
summary(mcmc0)
# model checking follows the course notes given in the MCMCglmm library
autocorr(mcmc0$Sol) # check the lowest non-zero lag values are lower than 0.1
autocorr(mcmc0$VCV) # bit high for residual variance (units)
# pdf("ModelComparisonRegions-plots.pdf", height=10, width=6)
plot(mcmc0$Sol) # summary plots for fixed effects
plot(mcmc0$VCV) # summary plots for random effects (and residual variance, "units")
posterior.mode(mcmc0$Sol) # get posterior distribution modes for each fixed effect
# I don't under stand how to do the nested continent/region structure in this function...

mcmc1 <- MCMCglmm(jack1.int ~ NPP_BP_mhw50_lg + Area_lg,  random=~GlobalBinCenter + Region, data=dats, family='poisson', nitt=30000)
summary(mcmc1)
autocorr(mcmc1$Sol)
autocorr(mcmc1$VCV)
plot(mcmc1$Sol)
plot(mcmc1$VCV)
mcmc1.1 <- MCMCglmm(jack1.int ~ NPP_BP_mhw50_lg + Area_lg + GlobalBinCenter, random=~GlobalBinCenter + Region, data=dats, family='poisson', nitt=30000)
summary(mcmc1.1)

# this next model does not work with default settings - need other priors... argh.
myprior<-list(R=list(V=1, nu=0), G=list(G1=list(V=diag(2), nu=2)))
mcmc2 <- MCMCglmm(jack1.int ~ NPP_BP_mhw50_lg + Area_lg,  random=~us(1+GlobalBinCenter):Region, data=dats, family='poisson', nitt=30000, prior=myprior)
summary(mcmc2)
autocorr(mcmc2$Sol)
autocorr(mcmc2$VCV)
plot(mcmc2$Sol)
plot(mcmc2$VCV)

mcmc4 <- MCMCglmm(jack1.int ~ NPP_BP_mhw50_lg + Area_lg + Time_lg + GlobalBinCenter, random=~Region, data=dats, family='poisson', nitt=30000)
summary(mcmc4)
autocorr(mcmc4$Sol)
autocorr(mcmc4$VCV)
plot(mcmc4$Sol)
plot(mcmc4$VCV)

mcmc5 <- MCMCglmm(jack1.int ~ NPP_BP_mhw50_lg + Area_lg + Time_lg,  random=~GlobalBinCenter + Region, data=dats, family='poisson', nitt=30000)
summary(mcmc5)
autocorr(mcmc5$Sol)
autocorr(mcmc5$VCV)
plot(mcmc5$Sol)
plot(mcmc5$VCV)

mcmc6 <- MCMCglmm(jack1.int ~ NPP_BP_mhw50_lg + Area_lg + Time_lg,  random=~us(1+GlobalBinCenter):Region, data=dats, family='poisson', nitt=30000, prior=myprior)
summary(mcmc6)
autocorr(mcmc6$Sol)
autocorr(mcmc6$VCV)
plot(mcmc6$Sol)
plot(mcmc6$VCV)
dev.off()

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

# final models for regions:
# run models and save residuals

# --   gamma diversity vs NPP from plants (best dates)
mm.bp <- glmer(jack1.int ~ NPP_BP_mhw50_lg + Area_lg + Time_lg + (GlobalBinCenter | Region), dats, family=poisson)
summary(mm.bp)
r.squaredGLMM(mm.bp)
dats$Mu.BP <- getME(mm.bp, "mu")
head(dats)
datsbp <- dats[,c(1,2,length(dats))]
datsbp$RegBin <- paste(datsbp$Region, datsbp$GlobalTimeBins, sep="_")

# --   gamma diversity vs NPP from plants (automatic dates)
dats <- savedats
dats <- dats[which(!is.na(dats$NPP_AP_mhw50_lg) & !is.na(dats$jack1)),] # dataset with complete cases
summary(dats)
mm.ap <- glmer(jack1.int ~ NPP_AP_mhw50_lg + Area_lg + Time_lg + (GlobalBinCenter | Region), dats, family=poisson)
summary(mm.ap)
r.squaredGLMM(mm.ap)
dats$Mu.AP <- getME(mm.ap, "mu")
datsap <- dats[,c(1,2,length(dats))]
datsap$RegBin <- paste(datsap$Region, datsap$GlobalTimeBins, sep="_")

# --   gamma diversity vs NPP from plants (best dates) with area accounting for tectonic extension of North America
dats <- savedats
dats <- dats[which(!is.na(dats$NPP_BP_mhw50_lg) & !is.na(dats$jack1)),] # dataset with complete cases
mm.bp2 <- glmer(jack1.int ~ NPP_BP_mhw50_lg + Area2_lg + Time_lg + (GlobalBinCenter | Continent), dats, family=poisson)
summary(mm.bp2)
r.squaredGLMM(mm.bp2)
datsbp2 <- dats[,c(1,2,length(dats))]
datsbp2$Mu.BP2 <- getME(mm.bp2, "mu")
datsbp2$RegBin <- paste(datsbp2$Region, datsbp2$GlobalTimeBins, sep="_")

# --   gamma diversity from cont average vs NPP from plants (best dates)
dats <- savedats
dats <- dats[which(!is.na(dats$NPP_BP_mhw50_lg) & !is.na(dats$jack1.conttime)),] # dataset with complete cases
mm.bp <- glmer(jack1.ct.int ~ NPP_BP_mhw50_lg + Area_lg + Time_lg + (GlobalBinCenter | Region), dats, family=poisson)
summary(mm.bp)
r.squaredGLMM(mm.bp)
dats$Mu.BP.ct <- getME(mm.bp, "mu")
datsbpct <- dats[,c(1,2,length(dats))]
datsbpct$RegBin <- paste(datsbpct$Region, datsbpct$GlobalTimeBins, sep="_")

# --   gamma diversity vs NPP from plants (automatic dates)
dats <- savedats
dats <- dats[which(!is.na(dats$NPP_AP_mhw50_lg) & !is.na(dats$jack1.conttime)),] # dataset with complete cases
summary(dats)
mm.ap <- glmer(jack1.ct.int ~ NPP_AP_mhw50_lg + Area_lg + Time_lg + (GlobalBinCenter | Region), dats, family=poisson)
summary(mm.ap)
dats$Mu.AP.ct <- getME(mm.ap, "mu")
datsapct <- dats[,c(1,2,length(dats))]
datsapct$RegBin <- paste(datsapct$Region, datsapct$GlobalTimeBins, sep="_")

# --   gamma diversity vs NPP from plants (best dates) with area accounting for tectonic extension of North America
dats <- savedats
dats <- dats[which(!is.na(dats$NPP_BP_mhw50_lg) & !is.na(dats$jack1.conttime)),] # dataset with complete cases
mm.bp2 <- glmer(jack1.ct.int ~ NPP_BP_mhw50_lg + Area2_lg + Time_lg + (GlobalBinCenter | Continent), dats, family=poisson)
summary(mm.bp2)
dats$Mu.BP2.ct <- getME(mm.bp2, "mu")
dats$RegBin <- paste(dats$Region, dats$GlobalTimeBins, sep="_")

# save the mu values
savedats$RegBin <- paste(savedats$Region, savedats$GlobalTimeBins, sep="_")
savedats$Mu.BP <- datsbp$Mu.BP[match(savedats$RegBin, datsbp$RegBin)]
savedats$Mu.BP2 <- datsbp2$Mu.BP2[match(savedats$RegBin, datsbp2$RegBin)]
savedats$Mu.AP <- datsap$Mu.AP[match(savedats$RegBin, datsap$RegBin)]
savedats$Mu.BP.ct <- datsbpct$Mu.BP.ct[match(savedats$RegBin, datsbpct$RegBin)]
savedats$Mu.BP2.ct <- dats$Mu.BP2.ct[match(savedats$RegBin, dats$RegBin)]
savedats$Mu.AP.ct <- datsapct$Mu.AP.ct[match(savedats$RegBin, datsapct$RegBin)]
savedats <- savedats[which(!is.na(savedats$GlobalTimeBins)),]
savedats[,24:length(savedats)]
names(savedats)
write.table(savedats, "Output_ModelDataRegions.txt", sep="\t", row.names=FALSE, col.names=TRUE)

