#### Regression on the global scale
#
# Global functions:
standardise <- function(x) (x-mean(x,na.rm = T))/sd(x, na.rm = T)
logit <- function(lo) log(lo/(1-lo)) 
invlogit <- function(x) 1/(1+exp(-x))
error_polygon <- function(ep,en,tstart,tend,tmid,color) {
  
  polygon( c(tstart, tmid, tend, tend, rev(tmid), tstart),
           c((ep)[1],ep, (ep)[length(ep)], (en)[length(en)], rev(en), (en)[1]),
           border = NA, col = color)
  
}
library(dplyr)
obiv <- readRDS("OBIS_bivalves_epifaunal_cleaned_2021_08_16_realms_and_grid.rds")

#  First prepare these tables (aggregate data on cell base)

############################################################################
#### Everything with Cells 6
gdat <- as.data.frame(matrix(NA,nrow=length(unique(obiv$cells6)), ncol = 1))
gdat <- 
  obiv %>%
  group_by(cells6) %>%
  dplyr::summarise(n = n(), nara = sum(isaragonite), lon = unique(gridlon6), lat = unique(gridlat6),
                   sstmean = mean(sstmean), sstrange = mean(sstrange),
                   sstmin = mean(sstmin), sstmax = mean(sstmax), realm = unique(realm),

                   n_ara_species = length(unique(species[which(isaragonite==1)])),n_species = length(unique(species))
  )

#table(gdat$realm)

evenness <- sapply(gdat$cells6, function(x) {
  spec = c(table(obiv$species[which(obiv$cells6 == x)]))
  specp = spec/sum(spec)
  -sum(specp*log(specp))/log(length(spec))
} )

gdat <- gdat %>% mutate(evenness = NA)
gdat$evenness <- evenness

dominance <- sapply(gdat$cells6, function(x) {
  spec = c(table(obiv$species[which(obiv$cells6 == x)]))
  specp = spec/sum(spec)
  max(specp)
} )
gdat <- gdat %>% mutate(dominance = NA)
gdat$dominance <- dominance


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

############################################################################
#### Everything with Cells 4
gdat4 <- as.data.frame(matrix(NA,nrow=length(unique(obiv$cells4)), ncol = 1))
gdat4 <- 
  obiv %>%
  group_by(cells4) %>%
  dplyr::summarise(n = n(), nara = sum(isaragonite), lon = unique(gridlon4), lat = unique(gridlat4),
                   sstmean = mean(sstmean), sstrange = mean(sstrange),
                   sstmin = mean(sstmin), sstmax = mean(sstmax), realm = unique(realm),

                   n_ara_species = length(unique(species[which(isaragonite==1)])),n_species = length(unique(species))
  )

#table(gdat4$realm)

evenness <- sapply(gdat4$cells4, function(x) {
  spec = c(table(obiv$species[which(obiv$cells4 == x)]))
  specp = spec/sum(spec)
  -sum(specp*log(specp))/log(length(spec))
} )

gdat4 <- gdat4 %>% mutate(evenness = NA)
gdat4$evenness <- evenness

dominance <- sapply(gdat4$cells4, function(x) {
  spec = c(table(obiv$species[which(obiv$cells4 == x)]))
  specp = spec/sum(spec)
  max(specp)
} )
gdat4 <- gdat4 %>% mutate(dominance = NA)
gdat4$dominance <- dominance

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

############################################################################
#### Everything with Cells 8
gdat8 <- as.data.frame(matrix(NA,nrow=length(unique(obiv$cells8)), ncol = 1))
gdat8 <- 
  obiv %>%
  group_by(cells8) %>%
  dplyr::summarise(n = n(), nara = sum(isaragonite), lon = unique(gridlon8), lat = unique(gridlat8),
                   sstmean = mean(sstmean), sstrange = mean(sstrange),
                   sstmin = mean(sstmin), sstmax = mean(sstmax), realm = unique(realm),
                   
                   n_ara_species = length(unique(species[which(isaragonite==1)])),n_species = length(unique(species))
  )

#table(gdat8$realm)

evenness <- sapply(gdat8$cells8, function(x) {
  spec = c(table(obiv$species[which(obiv$cells8 == x)]))
  specp = spec/sum(spec)
  -sum(specp*log(specp))/log(length(spec))
} )

gdat8 <- gdat8 %>% mutate(evenness = NA)
gdat8$evenness <- evenness

dominance <- sapply(gdat8$cells8, function(x) {
  spec = c(table(obiv$species[which(obiv$cells8 == x)]))
  specp = spec/sum(spec)
  max(specp)
} )
gdat8 <- gdat8 %>% mutate(dominance = NA)
gdat8$dominance <- dominance

############################################################################
####
####   Fig 2a
####
#png("Figs_06_2021/simple_reg_occ_g6_sstmean_stand.png", units = "in", width = 3.5, height = 3.5, res = 600)
windows(3.5,3.5)
#png("Revision_1/global_reg_occ_T_S_TS_tx.png", units = "in", width = 3.5, height = 3.5, res = 600)
#pdf("global_reg_occ_T_S_TS_tx.pdf", width = 3.5, height = 3.5)

dat2 <- gdat %>% filter(n_species > 1, evenness >= 0.5) #, nara > 0, nara/n < 1)


y <- cbind(dat2$nara,dat2$n-dat2$nara)

x1 <- standardise(dat2$sstmean)
x2 <- standardise(dat2$sstrange)

#x2 <- standardise(dat$sstrange)
# x1 <- standardise(datsamp$sstmin)
# x2 <- standardise(datsamp$sstmax)

glm1 <- (glm(y ~ x1*x2 , family = binomial(link = "logit")))
summary(glm1)


sizes <- (dat2$n)^0.11-0.45
par(mar = c(4,4,1,1), las = 1, mgp = c(2.5,0.75,0))
plot(x1,(dat2$nara/dat2$n), pch = 21, col = NA, bg = rgb(0,0,0,0.3), cex = sizes, xaxt = "n", xlab = expression("mean temperature ("*degree~"C)"), ylab = "P arag. occurrences")
axis(1, at = (c(0,10,20,30)-mean(dat2$sstmean))/sd(dat2$sstmean), labels = c(0,10,20,30) )
x1seq <- seq(min(x1),max(x1),0.1)


### Q 25

(quantile(dat2$sstrange,probs = 0.25)-mean(dat2$sstrange))/sd(dat2$sstrange)
newdata <- as.data.frame(cbind(x1seq,rep(-0.7767704 ,length(x1seq))))

colnames(newdata) = c("x1","x2")
predx1 <- predict.glm(glm1,newdata = newdata, se.fit = T)
predx1_05 <- invlogit(predx1$fit-2*predx1$se.fit)
predx1_95 <- invlogit(predx1$fit+2*predx1$se.fit)

error_polygon(predx1_05,predx1_95,x1seq[1],x1seq[length(x1seq)],x1seq,
             col = rgb(0,0.5,0.7,0.3))
points(x1seq,invlogit(predx1$fit),type = "l", lwd = 2, lty = 4,col = rgb(0,0.5,0.7,1))




### Q75

(quantile(dat2$sstrange, probs = 0.75)-mean(dat2$sstrange))/sd(dat2$sstrange)
newdata <- as.data.frame(cbind(x1seq,rep(0.6210254,length(x1seq))))

colnames(newdata) = c("x1","x2")
predx1 <- predict.glm(glm1,newdata = newdata, se.fit = T)
predx1_05 <- invlogit(predx1$fit-2*predx1$se.fit)
predx1_95 <- invlogit(predx1$fit+2*predx1$se.fit)

error_polygon(predx1_05,predx1_95,x1seq[1],x1seq[length(x1seq)],x1seq,
              col = rgb(0.7,0,0.5,0.3))
points(x1seq,invlogit(predx1$fit),type = "l", lwd = 2, lty = 2,col = rgb(0.7,0,0.5,1))




### Mean

newdata <- as.data.frame(cbind(x1seq,rep(0,length(x1seq))))
colnames(newdata) = c("x1","x2")
predx1 <- predict.glm(glm1,newdata = newdata, se.fit = T)
predx1_05 <- invlogit(predx1$fit-2*predx1$se.fit)
predx1_95 <- invlogit(predx1$fit+2*predx1$se.fit)

error_polygon(predx1_05,predx1_95,x1seq[1],x1seq[length(x1seq)],x1seq,
              col = rgb(0,0,0,0.3))
points(x1seq,invlogit(predx1$fit),type = "l", lwd = 2, lty = 1)

dev.off()



############################################################################
####
####   Fig 2b
####
#### T range on X

windows(3.5,3.5)
#png("Revision_1/global_reg_occ_T_S_TS_sx.png", units = "in", width = 3.5, height = 3.5, res = 600)
#pdf("global_reg_occ_T_S_TS_sx.pdf",  width = 3.5, height = 3.5)

dat2 <- gdat %>% filter(n_species > 1, evenness >= 0.5) #, nara > 0, nara/n < 1)


y <- cbind(dat2$nara,dat2$n-dat2$nara)

x1 <- standardise(dat2$sstmean)
x2 <- standardise(dat2$sstrange)

#x2 <- standardise(dat$sstrange)
# x1 <- standardise(datsamp$sstmin)
# x2 <- standardise(datsamp$sstmax)

glm1 <- (glm(y ~ x1*x2 , family = binomial(link = "logit")))
summary(glm1)


sizes <- (dat2$n)^0.11-0.45
#hist(sizes)
par(mar = c(4,4,1,1), las = 1, mgp = c(2.5,0.75,0))
plot(x2,(dat2$nara/dat2$n), pch = 21, col = NA, bg = rgb(0,0,0,0.3), cex = sizes, xlab = expression("seasonality ("*degree~"C)"), ylab = "P arag. occurrences",
     xaxs = "i", xlim = c((0-mean(dat2$sstrange))/sd(dat2$sstrange),(22-mean(dat2$sstrange))/sd(dat2$sstrange)), xaxt = "n")
axis(1, at = (c(0,5,10,15,20,30)-mean(dat2$sstrange))/sd(dat2$sstrange), labels = c(0,5,10,15,20,30) )
x1seq <- seq(min(x2),max(x2),0.1)


### Q 25

(quantile(dat2$sstmean,probs = 0.25)-mean(dat2$sstmean))/sd(dat2$sstmean)
newdata <- as.data.frame(cbind(rep(-0.7140352,length(x1seq)),x1seq))

colnames(newdata) = c("x1","x2")
predx1 <- predict.glm(glm1,newdata = newdata, se.fit = T)
predx1_05 <- invlogit(predx1$fit-2*predx1$se.fit)
predx1_95 <- invlogit(predx1$fit+2*predx1$se.fit)

error_polygon(predx1_05,predx1_95,x1seq[1],x1seq[length(x1seq)],x1seq,
              col = rgb(0,0.5,0.7,0.3))
points(x1seq,invlogit(predx1$fit),type = "l", lwd = 2, lty = 4,col = rgb(0,0.5,0.7,1))




### Q75

(quantile(dat2$sstmean, probs = 0.75)-mean(dat2$sstmean))/sd(dat2$sstmean)
newdata <- as.data.frame(cbind(rep(0.8099396,length(x1seq)),x1seq))

colnames(newdata) = c("x1","x2")
predx1 <- predict.glm(glm1,newdata = newdata, se.fit = T)
predx1_05 <- invlogit(predx1$fit-2*predx1$se.fit)
predx1_95 <- invlogit(predx1$fit+2*predx1$se.fit)

error_polygon(predx1_05,predx1_95,x1seq[1],x1seq[length(x1seq)],x1seq,
              col = rgb(0.7,0,0.5,0.3))
points(x1seq,invlogit(predx1$fit),type = "l", lwd = 2, lty = 2,col = rgb(0.7,0,0.5,1))

a <- min(x1)
b <- min(x2)
invlogit(-0.99073+0.41283 *a  + 0.01656    * b + 0.29746    *a*b )


### Mean

newdata <- as.data.frame(cbind(rep(0,length(x1seq)),x1seq))
colnames(newdata) = c("x1","x2")
predx1 <- predict.glm(glm1,newdata = newdata, se.fit = T)
predx1_05 <- invlogit(predx1$fit-2*predx1$se.fit)
predx1_95 <- invlogit(predx1$fit+2*predx1$se.fit)

error_polygon(predx1_05,predx1_95,x1seq[1],x1seq[length(x1seq)],x1seq,
              col = rgb(0,0,0,0.3))
points(x1seq,invlogit(predx1$fit),type = "l", lwd = 2, lty = 1)


dev.off()



############################################################################
####
####   Fig 2c
####
#### Figure with standardised global regression coefficients
#
#
windows(4.3,3.2)

#png("Revision_1/global_reg_occ_coefficients.png", units = "in", width = 4.3, height = 3.2, res = 600)
pdf("global_reg_occ_coefficients.pdf",  width = 4.3, height = 3.2)

baseat <- c(1,2,3,3.85,4.7)
dat <- gdat %>% filter(n_species > 1, evenness >= 0.5)



#y <- cbind(dat$n_ara_species,dat$n_species)

y <- cbind(dat$nara,dat$n-dat$nara)

x1 <- standardise(dat$sstmean)
x2 <- standardise(dat$sstrange)
#x1 <- standardise(dat$sstmin)
#x2 <- standardise(dat$sstmax)

glm1 <- (glm(y ~ x1 , family = binomial(link = "logit")))
glm2 <- (glm(y ~ x2 , family = binomial(link = "logit")))
glm12 <- (glm(y ~ x1 * x2 , family = binomial(link = "logit")))

s1 <- summary(glm1)
s2 <- summary(glm2)
s12 <- summary(glm12)


mod1 <- -c(-0.3,-0.15,0,0.15,0.3 )
par(mar = c(4,4,1,1), las = 1, mgp = c(2.5,0.75,0))
plot(0,0,type= "n", xlim = c(0.5,5.2), xaxt = "n", xlab = NA, ylab = "slope", xaxs = "i", yaxs = "i", ylim = c(-0.58,0.74), yaxt = "n")
axis(2,c(-1,-0.2,0,0.2,0.4,0.6,1))


legend("bottomleft", pch = c(1,21,21,25,22), pt.cex = c(1.45,1.1,1.1,1,1.1), legend = c(expression("\xd8"~"= 247 km, "~italic("R")>=" 2"),  # "~italic("J'")~">= 0.5"),
                                                                                        expression("\xd8"~"= 247 km, "~italic("R")>=" 10"),
                                                                                        expression("\xd8"~"= 247 km, "~italic("R")>=" 2, "~italic("D")>=" 0.75 "), #"~italic("J'")~">= 0.5"),
                                                                                        expression("\xd8"~"= 555 km, "~italic("R")>=" 2"), #  "~italic("J'")~">= 0.5"),
                                                                                        expression("\xd8"~"= 139 km, "~italic("R")>=" 2")), #  "~italic("J'")~">= 0.5")),
       cex = 0.71, box.col = NA, bty = "o", bg = "white", pt.bg = rgb(0,0,0,0.5), col = c("black", NA, "black",NA, NA),ncol = 2)

points(c(1.5,1.5), c(-0.25,0.73), type = "l")
points(c(2.5,2.5), c(-0.25,0.73), type = "l")
points(c(3.425,3.425), c(-0.25,0.73), type = "l", col = "grey75")
points(c(4.275,4.275), c(-0.25,0.73), type = "l", col = "grey75")

abline(h = c(-0.25,-0.58), lwd = c(1,1.5))


abline(h = 0, lty = 3)
xat <- baseat-mod1[1]
col1 <- rgb(0,0,0,0.5)
cex1 <- 1.45
pch1 <- 1
points(xat[1],s1$coefficients[2,1], pch = pch1, cex = cex1, col = "black", bg = col1)
points(rep(xat[1],2), s1$coefficients[2,1]+c(-2,2)*s1$coefficients[2,2], type = "l", lwd = 2)

ind <- 2
points(xat[ind],s2$coefficients[2,1], pch = pch1, cex = cex1, col = "black", bg = col1)
points(rep(xat[ind],2), s2$coefficients[2,1]+c(-2,2)*s2$coefficients[2,2], type = "l", lwd = 2)

ind <- 3
points(xat[ind],s12$coefficients[2,1], pch = pch1, cex = cex1, col = "black", bg = col1)
points(rep(xat[ind],2), s12$coefficients[2,1]+c(-2,2)*s12$coefficients[2,2], type = "l", lwd = 2)

ind <- 4
points(xat[ind],s12$coefficients[3,1], pch = pch1, cex = cex1, col = "black", bg = col1)
points(rep(xat[ind],2), s12$coefficients[3,1]+c(-2,2)*s12$coefficients[3,2], type = "l", lwd = 2)

ind <- 5
points(xat[ind],s12$coefficients[4,1], pch = pch1, cex = cex1, col = "black", bg = col1)
points(rep(xat[ind],2), s12$coefficients[4,1]+c(-2,2)*s12$coefficients[4,2], type = "l", lwd = 2)

###

dat <- gdat %>% filter(n_species > 9, evenness >= 0.5)


#y <- cbind(dat$n_ara_species,dat$n_species)

y <- cbind(dat$nara,dat$n-dat$nara)

x1 <- standardise(dat$sstmean)
x2 <- standardise(dat$sstrange)
#x1 <- standardise(dat$sstmin)
#x2 <- standardise(dat$sstmax)

glm1 <- (glm(y ~ x1 , family = binomial(link = "logit")))
glm2 <- (glm(y ~ x2 , family = binomial(link = "logit")))
glm12 <- (glm(y ~ x1 * x2 , family = binomial(link = "logit")))

s1 <- summary(glm1)
s2 <- summary(glm2)
s12 <- summary(glm12)

xat <- baseat-mod1[2]
col1 <- rgb(0,0,0,0.5)
cex1 <- 1.1
pch1 <- 21
points(xat[1],s1$coefficients[2,1], pch = pch1, cex = cex1, col = NA, bg = col1)
points(rep(xat[1],2), s1$coefficients[2,1]+c(-2,2)*s1$coefficients[2,2], type = "l", lwd = 2)

ind <- 2
points(xat[ind],s2$coefficients[2,1], pch = pch1, cex = cex1, col = NA, bg = col1)
points(rep(xat[ind],2), s2$coefficients[2,1]+c(-2,2)*s2$coefficients[2,2], type = "l", lwd = 2)

ind <- 3
points(xat[ind],s12$coefficients[2,1], pch = pch1, cex = cex1, col = NA, bg = col1)
points(rep(xat[ind],2), s12$coefficients[2,1]+c(-2,2)*s12$coefficients[2,2], type = "l", lwd = 2)

ind <- 4
points(xat[ind],s12$coefficients[3,1], pch = pch1, cex = cex1, col = NA, bg = col1)
points(rep(xat[ind],2), s12$coefficients[3,1]+c(-2,2)*s12$coefficients[3,2], type = "l", lwd = 2)

ind <- 5
points(xat[ind],s12$coefficients[4,1], pch = pch1, cex = cex1, col = NA, bg = col1)
points(rep(xat[ind],2), s12$coefficients[4,1]+c(-2,2)*s12$coefficients[4,2], type = "l", lwd = 2)


###


dat <- gdat %>% filter(n_species > 1, dominance <= 0.75)


#y <- cbind(dat$n_ara_species,dat$n_species)

y <- cbind(dat$nara,dat$n-dat$nara)

x1 <- standardise(dat$sstmean)
x2 <- standardise(dat$sstrange)
#x1 <- standardise(dat$sstmin)
#x2 <- standardise(dat$sstmax)

glm1 <- (glm(y ~ x1 , family = binomial(link = "logit")))
glm2 <- (glm(y ~ x2 , family = binomial(link = "logit")))
glm12 <- (glm(y ~ x1 * x2 , family = binomial(link = "logit")))

s1 <- summary(glm1)
s2 <- summary(glm2)
s12 <- summary(glm12)

xat <- baseat-mod1[3]
col1 <- rgb(0,0,0,0.5)
cex1 <- 1.1
pch1 <- 21
points(xat[1],s1$coefficients[2,1], pch = pch1, cex = cex1, col = "black", bg = col1)
points(rep(xat[1],2), s1$coefficients[2,1]+c(-2,2)*s1$coefficients[2,2], type = "l", lwd = 2)

ind <- 2
points(xat[ind],s2$coefficients[2,1], pch = pch1, cex = cex1, col = "black", bg = col1)
points(rep(xat[ind],2), s2$coefficients[2,1]+c(-2,2)*s2$coefficients[2,2], type = "l", lwd = 2)

ind <- 3
points(xat[ind],s12$coefficients[2,1], pch = pch1, cex = cex1, col = "black", bg = col1)
points(rep(xat[ind],2), s12$coefficients[2,1]+c(-2,2)*s12$coefficients[2,2], type = "l", lwd = 2)

ind <- 4
points(xat[ind],s12$coefficients[3,1], pch = pch1, cex = cex1, col = "black", bg = col1)
points(rep(xat[ind],2), s12$coefficients[3,1]+c(-2,2)*s12$coefficients[3,2], type = "l", lwd = 2)

ind <- 5
points(xat[ind],s12$coefficients[4,1], pch = pch1, cex = cex1, col = "black", bg = col1)
points(rep(xat[ind],2), s12$coefficients[4,1]+c(-2,2)*s12$coefficients[4,2], type = "l", lwd = 2)


###



dat <- gdat4 %>% filter(n_species > 1, evenness >= 0.5)


#y <- cbind(dat$n_ara_species,dat$n_species)

y <- cbind(dat$nara,dat$n-dat$nara)

x1 <- standardise(dat$sstmean)
x2 <- standardise(dat$sstrange)
#x1 <- standardise(dat$sstmin)
#x2 <- standardise(dat$sstmax)

glm1 <- (glm(y ~ x1 , family = binomial(link = "logit")))
glm2 <- (glm(y ~ x2 , family = binomial(link = "logit")))
glm12 <- (glm(y ~ x1 * x2 , family = binomial(link = "logit")))

s1 <- summary(glm1)
s2 <- summary(glm2)
s12 <- summary(glm12)

xat <- baseat-mod1[4]
col1 <- rgb(0,0,0,0.5)
cex1 <- 1
pch1 <- 25
points(xat[1],s1$coefficients[2,1], pch = pch1, cex = cex1, col = NA, bg = col1)
points(rep(xat[1],2), s1$coefficients[2,1]+c(-2,2)*s1$coefficients[2,2], type = "l", lwd = 2)

ind <- 2
points(xat[ind],s2$coefficients[2,1], pch = pch1, cex = cex1, col = NA, bg = col1)
points(rep(xat[ind],2), s2$coefficients[2,1]+c(-2,2)*s2$coefficients[2,2], type = "l", lwd = 2)

ind <- 3
points(xat[ind],s12$coefficients[2,1], pch = pch1, cex = cex1, col = NA, bg = col1)
points(rep(xat[ind],2), s12$coefficients[2,1]+c(-2,2)*s12$coefficients[2,2], type = "l", lwd = 2)

ind <- 4
points(xat[ind],s12$coefficients[3,1], pch = pch1, cex = cex1, col = NA, bg = col1)
points(rep(xat[ind],2), s12$coefficients[3,1]+c(-2,2)*s12$coefficients[3,2], type = "l", lwd = 2)

ind <- 5
points(xat[ind],s12$coefficients[4,1], pch = pch1, cex = cex1, col = NA, bg = col1)
points(rep(xat[ind],2), s12$coefficients[4,1]+c(-2,2)*s12$coefficients[4,2], type = "l", lwd = 2)

###

dat <- gdat8 %>% filter(n_species > 1, evenness >= 0.5)


#y <- cbind(dat$n_ara_species,dat$n_species)

y <- cbind(dat$nara,dat$n-dat$nara)

x1 <- standardise(dat$sstmean)
x2 <- standardise(dat$sstrange)
#x1 <- standardise(dat$sstmin)
#x2 <- standardise(dat$sstmax)

glm1 <- (glm(y ~ x1 , family = binomial(link = "logit")))
glm2 <- (glm(y ~ x2 , family = binomial(link = "logit")))
glm12 <- (glm(y ~ x1 * x2 , family = binomial(link = "logit")))

s1 <- summary(glm1)
s2 <- summary(glm2)
s12 <- summary(glm12)

xat <- baseat-mod1[5]
col1 <- rgb(0,0,0,0.5)
cex1 <- 1.1
pch1 <- 22
points(xat[1],s1$coefficients[2,1], pch = pch1, cex = cex1, col = NA, bg = col1)
points(rep(xat[1],2), s1$coefficients[2,1]+c(-2,2)*s1$coefficients[2,2], type = "l", lwd = 2)

ind <- 2
points(xat[ind],s2$coefficients[2,1], pch = pch1, cex = cex1, col = NA, bg = col1)
points(rep(xat[ind],2), s2$coefficients[2,1]+c(-2,2)*s2$coefficients[2,2], type = "l", lwd = 2)

ind <- 3
points(xat[ind],s12$coefficients[2,1], pch = pch1, cex = cex1, col = NA, bg = col1)
points(rep(xat[ind],2), s12$coefficients[2,1]+c(-2,2)*s12$coefficients[2,2], type = "l", lwd = 2)

ind <- 4
points(xat[ind],s12$coefficients[3,1], pch = pch1, cex = cex1, col = NA, bg = col1)
points(rep(xat[ind],2), s12$coefficients[3,1]+c(-2,2)*s12$coefficients[3,2], type = "l", lwd = 2)

ind <- 5
points(xat[ind],s12$coefficients[4,1], pch = pch1, cex = cex1, col = NA, bg = col1)
points(rep(xat[ind],2), s12$coefficients[4,1]+c(-2,2)*s12$coefficients[4,2], type = "l", lwd = 2)

axis(1,at = baseat, labels = expression(italic("T"), italic("S"), italic("T"), italic("S"), italic("T")~" \u00D7"~italic("S")))



abline(v = c(0.5,5.2), lwd = c(1.5,1.5))

dev.off()
############################################################################


############################################################################
####
####   Get D-squared for the regressions
####   https://modtools.wordpress.com/2013/08/14/dsquared/
#### 

Dsquared <- function(model = NULL, 
                     obs = NULL, 
                     pred = NULL, 
                     family = NULL, # needed only when 'model' not provided
                     adjust = FALSE, 
                     npar = NULL) { # needed only when 'model' not provided
  # version 1.4 (31 Aug 2015)
  
  model.provided <- ifelse(is.null(model), FALSE, TRUE)
  
  if (model.provided) {
    if (!("glm" %in% class(model))) stop ("'model' must be of class 'glm'.")
    if (!is.null(pred)) message("Argument 'pred' ignored in favour of 'model'.")
    if (!is.null(obs)) message("Argument 'obs' ignored in favour of 'model'.")
    obs <- model$y
    pred <- model$fitted.values
    
  } else { # if model not provided
    if (is.null(obs) | is.null(pred)) stop ("You must provide either 'obs' and 'pred', or a 'model' object of class 'glm'.")
    if (length(obs) != length(pred)) stop ("'obs' and 'pred' must be of the same length (and in the same order).")
    if (is.null(family)) stop ("With 'obs' and 'pred' arguments (rather than a model object), you must also specify one of two model family options: 'binomial' or 'poisson' (in quotes).")
    else if (!is.character(family)) stop ("Argument 'family' must be provided as character (i.e. in quotes: 'binomial' or 'poisson').")
    else if (length(family) != 1 | !(family %in% c("binomial", "poisson"))) stop ("'family' must be either 'binomial' or 'poisson' (in quotes).")
    
    if (family == "binomial") {
      if (any(!(obs %in% c(0, 1)) | pred < 0 | pred > 1)) stop ("'binomial' family implies that 'obs' data should be binary (with values 0 or 1) and 'pred' data should be bounded between 0 and 1.")
      link <- log(pred / (1 - pred))  # logit
    }  # end if binomial
    
    else if (family == "poisson") {
      if (any(obs %%1 != 0)) stop ("'poisson' family implies that 'obs' data should consist of whole numbers.")
      link <- log(pred)
    }  # end if poisson
    
    model <- glm(obs ~ link, family = family)
  }  # end if model not provided
  
  D2 <- (model$null.deviance - model$deviance) / model$null.deviance
  
  if (adjust) {
    if (model.provided) {
      n <- length(model$y)
      #p <- length(model$coefficients)
      p <- attributes(logLik(model))$df
    } else {
      if (is.null(npar)) stop ("Adjusted D-squared from 'obs' and 'pred' values (rather than a model object) requires specifying the number of parameters in the underlying model ('npar').")
      n <- length(na.omit(obs))
      p <- npar
    }  # end if model.provided else
    
    D2 <- 1 - ((n - 1) / (n - p)) * (1 - D2)
  }  # end if adjust
  
  return (D2)
}




dat <- gdat %>% filter(n_species > 1, evenness >= 0.5)


#y <- cbind(dat$n_ara_species,dat$n_species)

y <- cbind(dat$nara,dat$n-dat$nara)

x1 <- standardise(dat$sstmean)
x2 <- standardise(dat$sstrange)
#x1 <- standardise(dat$sstmin)
#x2 <- standardise(dat$sstmax)

glm1 <- (glm(y ~ x1 , family = binomial(link = "logit")))
glm2 <- (glm(y ~ x2 , family = binomial(link = "logit")))
glm12 <- (glm(y ~ x1 * x2 , family = binomial(link = "logit")))

Dsquared(glm1, adjust = FALSE)

