#### Regression on the scale of biogeographic realms
#
obiv <- readRDS("OBIS_bivalves_epifaunal_cleaned_2021_08_16_realms_and_grid.rds")
library(dplyr)
library(R2jags)
library(ggmcmc)

# 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)
  
}
#
#### Prepare data
#
### Realm colours
### Realm colours:
alph <- 1*255
realmcolours1 <- (c(rgb(201,87,109,alph, maxColorValue = 255),
                    rgb(225,148,66,alph, maxColorValue = 255),
                    
                    rgb(129,94,205,alph, maxColorValue = 255),
                    rgb(174,173,67,alph, maxColorValue = 255),
                    rgb(81,180,141,alph, maxColorValue = 255),
                    
                    rgb(81,127,56,alph, maxColorValue = 255),
                    rgb(75,175,208,alph, maxColorValue = 255),
                    rgb(101,182,75,alph, maxColorValue = 255),
                    
                    rgb(200,101,185,alph, maxColorValue = 255),
                    
                    rgb(93,129,216,alph, maxColorValue = 255),
                    rgb(153,114,56,alph, maxColorValue = 255),
                    rgb(205,81,62,alph, maxColorValue = 255))[12:1])[c(1,7,6,3,11,10,
                                                                       12,2,8,9,5,4)]

alph <- 0.35*255

realmcolours2 <- (c(rgb(201,87,109,alph, maxColorValue = 255),
                    rgb(225,148,66,alph, maxColorValue = 255),
                    
                    rgb(129,94,205,alph, maxColorValue = 255),
                    rgb(174,173,67,alph, maxColorValue = 255),
                    rgb(81,180,141,alph, maxColorValue = 255),
                    
                    rgb(81,127,56,alph, maxColorValue = 255),
                    rgb(75,175,208,alph, maxColorValue = 255),
                    rgb(101,182,75,alph, maxColorValue = 255),
                    
                    rgb(200,101,185,alph, maxColorValue = 255),
                    
                    rgb(93,129,216,alph, maxColorValue = 255),
                    rgb(153,114,56,alph, maxColorValue = 255),
                    rgb(205,81,62,alph, maxColorValue = 255))[12:1])


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

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



model_interaction <- function() {
  
  for(j in 1:Nrealm) {
    for(i in 1:Nperrealm[j]) {
      r[i,j] ~ dbin(p[i,j],n[i,j])
      logit(p[i,j])  <- alpha[j] + beta1[j]*x1[i,j] + beta2[j]*x2[i,j] + beta3[j]*x1[i,j]*x2[i,j]
    }
  }
  
  for(j in 1:Nrealm)  {
    alpha[j] ~ dnorm(alpha_mu,alpha_tau)
    beta1[j] ~ dnorm(beta1_mu,beta1_tau)
    beta2[j] ~ dnorm(beta2_mu,beta2_tau)
    beta3[j] ~ dnorm(beta3_mu,beta3_tau)
    
  }
  alpha_mu ~ dnorm(0,0.01)
  beta1_mu ~ dnorm(0,0.01)
  beta2_mu ~ dnorm(0,0.01)
  beta3_mu ~ dnorm(0,0.01)
  
  alpha_log_tau ~ dnorm(0,0.01)
  beta1_log_tau ~ dnorm(0,0.01)
  beta2_log_tau ~ dnorm(0,0.01)
  beta3_log_tau ~ dnorm(0,0.01)
  
  alpha_tau <- exp(alpha_log_tau) 
  beta1_tau <- exp(beta1_log_tau) 
  beta2_tau <- exp(beta2_log_tau) 
  beta3_tau <- exp(beta3_log_tau) 
  
}

# Filter to include only the eight realms with decent sampling
dat <- gdat %>% filter(n_species > 1, evenness >= 0.5, realm %in% c(2,3,4,6,8,9,11,12))
#dat <- gdat %>% filter(n_species > 1, evenness >= 0.5, realm %in% c(3,6,8,11))

table(dat$realm)

Nrealm <- 8
#N_10  <- c(3,  6,  8, 11)
N_10 <- c(2,3,4,6,8,9,11,12)
#N_10 <- c(3,6,8,11)

#dat <- gdat %>% filter(n_species > 4, evenness >= 0.5, realm %in% N_10)

Nperrealm <- c(table(dat$realm))


x1 <- matrix(NA,nrow = max(Nperrealm), ncol = Nrealm)
for(i in 1:Nrealm) x1[1:Nperrealm[i],i] <- standardise(dat$sstmean[which(dat$realm==N_10[i])])

# x1[1,5] <- 0
x2 <- matrix(NA,nrow = max(Nperrealm), ncol = Nrealm)
for(i in 1:Nrealm) x2[1:Nperrealm[i],i] <- standardise(dat$sstrange[which(dat$realm==N_10[i])])

r <- matrix(NA,nrow = max(Nperrealm), ncol = Nrealm)
for(i in 1:Nrealm) r[1:Nperrealm[i],i] <- dat$nara[which(dat$realm==N_10[i])]

n <- matrix(NA,nrow = max(Nperrealm), ncol = Nrealm)
for(i in 1:Nrealm) n[1:Nperrealm[i],i] <- dat$n[which(dat$realm==N_10[i])]

data_jags <- list(x1 = x1, x2 = x2, r = r, n = n,
                  Nrealm = Nrealm, Nperrealm = Nperrealm)

### Run analysis
reg_t_s_ts_8r <- jags.parallel(data = data_jags, 
                                                       parameters.to.save = c("alpha",
                                                                              "beta1",
                                                                              "beta2",
                                                                              "beta3",
                                                                              "alpha_mu",
                                                                              "alpha_tau",
                                                                              "beta1_mu",
                                                                              "beta1_tau",
                                                                              "beta2_mu",
                                                                              "beta2_tau",
                                                                              "beta3_mu",
                                                                              "beta3_tau"
                                                       ),
                                                       n.iter = 15000,
                                                       n.burnin = 5000,
                                                       n.thin = 1,
                                                       n.chains =  3, # Other values set at default (for simplicity)
                                                       model.file = model_interaction)

### Check the posterior distribution and Markov chains for individual parameters: 
reg_t_s_ts_8r.mcmcm <- as.mcmc(reg_t_s_ts_8r)
reg_t_s_ts_8rg <- ggs(reg_t_s_ts_8r.mcmcm)
ggs_traceplot(reg_t_s_ts_8rg, family = "beta3_mu")
ggs_density(reg_t_s_ts_8rg, family = "beta3_mu")

# Do analogous regressions with just T and just R
#
# For the plots, we will simply load the results of those analyses

#########################################################
###
### Fig 3d



reg2 <- readRDS("reg_t_s_ts_8r.rds")
reg1 <- readRDS("reg_t_8r.rds")
rp1 <- reg1$BUGSoutput$sims.list
rp2 <- reg2$BUGSoutput$sims.list



N_10 <- c(2,3,4,6,8,9,11,12)


### X1 X2 reg

windows(5.2,3.1)
#png("Revision_1/hierarch_reg_occ_g6_12_slopes_stand_meanT_8r.png", units = "in", width = 5.2, height = 3.1, res = 600)
#pdf("hierarch_reg_occ_g6_12_slopes_stand_meanT_8r.pdf", width = 5.2, height = 3.1)

par(mar=c(4.1,4.1,1,1), las = 1, mgp = c(2.4,0.75,0), mfrow = c(1,1))
plot(0,0,type = "n", xlim = c(-0.5,12.5), ylim = c(-1,3), xaxs = "i", yaxs = "i",
     xlab = "Realm", ylab = "temperature slope", xaxt = "n")

axis(1,at = c(-1,0,N_10,13), labels = c(NA,expression(mu[beta]),N_10,NA))
abline(h = 0)

#abline(h = mean(rp1$beta1_mu),col = rgb(0.3,0.3,0.3,1), lwd = 1.5, lty = 2)
points(c(0,0)-.2,c(quantile(rp1$beta1_mu, probs = c(0.025,0.975))), col = rgb(0.3,0.3,0.3,1), lwd = 2.4, type = "l")
points(c(0)-.2,mean(rp1$beta1_mu), col = rgb(0.3,0.3,0.3,1), pch = 21, cex = 1.5, bg = rgb(0.3,0.3,0.3,0.4))

for(i in 1:8) {
  region <- N_10[i]
  points(c(region)-.2,mean(rp1$beta1[,i]), col = realmcolours1[region], pch = 21, cex = 1.2, lwd = 1.2, bg = realmcolours2[region])
  points(c(region,region)-.2,c(quantile(rp1$beta1[,i],probs = c(0.025)),quantile(rp1$beta1[,i],probs = c(0.975))), type = "l", lwd = 2, col = realmcolours1[region])
  
}

points(c(0,0)+.2,c(quantile(rp2$beta1_mu, probs = c(0.025,0.975))), col = rgb(0.3,0.3,0.3,1), lwd = 2.4, type = "l")
points(c(0)+.2,mean(rp2$beta1_mu), col = rgb(0.3,0.3,0.3,1), pch = 22, cex = 1.4, bg = rgb(0.3,0.3,0.3,0.4))

for(i in 1:8) {
 region <- N_10[i]
    points(c(region)+.2,mean(rp2$beta1[,i]), col = realmcolours1[region], pch = 22, cex = 1.1, lwd = 1.2, bg = realmcolours2[region])
  points(c(region,region)+.2,c(quantile(rp2$beta1[,i],probs = c(0.025)),quantile(rp2$beta1[,i],probs = c(0.975))), type = "l", lwd = 2, col = realmcolours1[region])
  
}

dev.off()



#########################################################
###
### Fig 3e

reg2 <- readRDS("reg_t_s_ts_8r.rds")
reg1 <- readRDS("reg_s_8r.rds")
rp1 <- reg1$BUGSoutput$sims.list
rp2 <- reg2$BUGSoutput$sims.list

### X1 X2 reg

 windows(5.2,3.1)
#cairo_ps("Regression_biv_occ_min_Temperature_x.eps", width = 3.5, height  = 3.5, fallback_resolution = 600)
#png("Revision_1/hierarch_reg_occ_g6_12_slopes_stand_Trange_8r.png", units = "in", width = 5.2, height = 3.1, res = 600)
 #pdf("hierarch_reg_occ_g6_12_slopes_stand_Trange_8r.pdf", width = 5.2, height = 3.1)
 
par(mar=c(4.1,4.1,1,1), las = 1, mgp = c(2.4,0.75,0), mfrow = c(1,1))
plot(0,0,type = "n", xlim = c(-0.5,12.5), ylim = c(-1.1,1.6), xaxs = "i", yaxs = "i",
     xlab = "Realm", ylab = "seasonality slope", xaxt = "n")

axis(1,at = seq(-1,13,1), labels = c(NA,expression(mu[beta]),1:12,NA))
abline(h = 0)

#abline(h = mean(rp1$beta1_mu),col = rgb(0.3,0.3,0.3,1), lwd = 1.5, lty = 2)

points(c(0,0)-.2,c(quantile(rp1$beta1_mu, probs = c(0.025,0.975))), col = rgb(0.3,0.3,0.3,1), lwd = 2.4, type = "l")
points(c(0)-.2,mean(rp1$beta1_mu), col = rgb(0.3,0.3,0.3,1), pch = 21, cex = 1.5, bg = rgb(0.3,0.3,0.3,0.4))

for(i in 1:8) {
  region <- N_10[i]
  points(c(region)-.2,mean(rp1$beta1[,i]), col = realmcolours1[region], pch = 21, cex = 1.2, lwd = 1.2, bg = realmcolours2[region])
  points(c(region,region)-.2,c(quantile(rp1$beta1[,i],probs = c(0.025)),quantile(rp1$beta1[,i],probs = c(0.975))), type = "l", lwd = 2, col = realmcolours1[region])
  
}

points(c(0,0)+.2,c(quantile(rp2$beta2_mu, probs = c(0.025,0.975))), col = rgb(0.3,0.3,0.3,1), lwd = 2.4, type = "l")
points(c(0)+.2,mean(rp2$beta2_mu), col = rgb(0.3,0.3,0.3,1), pch = 22, cex = 1.4, bg = rgb(0.3,0.3,0.3,0.4))

for(i in 1:8) {
  region <- N_10[i]
  points(c(region)+.2,mean(rp2$beta2[,i]), col = realmcolours1[region], pch = 22, cex = 1.1, lwd = 1.2, bg = realmcolours2[region])
  points(c(region,region)+.2,c(quantile(rp2$beta2[,i],probs = c(0.025)),quantile(rp2$beta2[,i],probs = c(0.975))), type = "l", lwd = 2, col = realmcolours1[region])
  
}

dev.off()

#########################################################
###
### Fig 3f

windows(5.2,3.1)
#cairo_ps("Regression_biv_occ_min_Temperature_x.eps", width = 3.5, height  = 3.5, fallback_resolution = 600)
#png("Revision_1/hierarch_reg_occ_g6_12_slopes_stand_interaction_8r.png", units = "in", width = 5.2, height = 3.1, res = 600)
#pdf("hierarch_reg_occ_g6_12_slopes_stand_interaction_8r.pdf", width = 5.2, height = 3.1)

par(mar=c(4.1,4.1,1,1), las = 1, mgp = c(2.4,0.75,0), mfrow = c(1,1))
plot(0,0,type = "n", xlim = c(-0.5,12.5), ylim = c(-1,1), xaxs = "i", yaxs = "i",
     xlab = "Realm", ylab = expression(italic("T")~" \u00D7"~italic("S")~"slope"), xaxt = "n")

axis(1,at = seq(-1,13,1), labels = c(NA,expression(mu[beta]),1:12,NA))
abline(h = 0)


points(c(0,0),c(quantile(rp2$beta3_mu, probs = c(0.025,0.975))), col = rgb(0.3,0.3,0.3,1), lwd = 2.4, type = "l")
points(c(0),mean(rp2$beta3_mu), col = rgb(0.3,0.3,0.3,1), pch = 22, cex = 1.4, bg = rgb(0.3,0.3,0.3,0.4))

offset1 <- rep(c(0,1),4)
for(i in 1:8) {
  region <- N_10[i]
  points(c(region),mean(rp2$beta3[,i]), col = realmcolours1[region], pch = 22, cex = 1.1, lwd = 1.2, bg = realmcolours2[region])
  points(c(region,region),c(quantile(rp2$beta3[,i],probs = c(0.025)),quantile(rp2$beta3[,i],probs = c(0.975))), type = "l", lwd = 2, col = realmcolours1[region])
  text(region,0.9+offset1[i]*-0.2,labels = sprintf("%.2f",round((cor.test(x1[,i],x2[,i]))$estimate,2)), cex = 0.825)
}

text(0.75,0.8,"r = ", cex = 0.9)

dev.off()


#########################################################
###
### Fig 3b


#### Visualise the overall hierarchical regression

rp2 <- reg2$BUGSoutput$sims.list

#### T  on X

#### T  on X

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

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


x1 <- matrix(NA,nrow = max(Nperrealm), ncol = Nrealm)
for(i in 1:Nrealm) x1[1:Nperrealm[i],i] <- standardise(dat$sstmean[which(dat$realm==N_10[i])])

# x1[1,5] <- 0
x2 <- matrix(NA,nrow = max(Nperrealm), ncol = Nrealm)
for(i in 1:Nrealm) x2[1:Nperrealm[i],i] <- standardise(dat$sstrange[which(dat$realm==N_10[i])])

x1 <- na.omit(c(x1))
x2 <- na.omit(c(x2))


sizes <- (dat2$n)^0.11-0.45
par(mar = c(4,4,1,1), las = 1, mgp = c(2.5,0.75,0))
plot(0,0, pch = 21, col = NA, bg = rgb(0,0,0,0.3), cex = sizes, xlab = expression("standardised temperature"), ylab = "P arag. occurrences",
     xaxs = "i", xlim = c(min(x1)-0.1,max(x1)+0.1), type = "n", ylim = c(0,1))
#axis(1, at = (c(-10,0,10,20,30,40)-mean(dat2$sstmean))/sd(dat2$sstmean), labels = c(-10,0,10,20,30,40) )
x1seq <- seq(min(x1),max(x1),0.1)



### Mean
bigN <- nrow(rp2$alpha_mu)
bigN_index <- seq(1,bigN,1)
bigN <- length(bigN_index)
pred1 <- matrix(NA,nrow = bigN,ncol = length(x1seq))
for(i in 1:bigN)  pred1[i,] <- invlogit((rp2$alpha_mu[i]) + rp2$beta1_mu[i]*x1seq + rp2$beta2_mu[i]*0 + rp2$beta3_mu[i]*0*x1seq)

low1 <- apply(pred1,2,function(x) quantile(x,0.025))
high1 <- apply(pred1,2,function(x) quantile(x,0.975))
med1 <- apply(pred1,2,function(x) quantile(x,0.5))
mean1 <- apply(pred1,2,function(x) mean(x))

error_polygon(low1,high1,x1seq[1],x1seq[length(x1seq)],x1seq,
              col = rgb(0,0,0,0.25))
points(x1seq,invlogit(mean(rp2$alpha_mu)+mean(rp2$beta1_mu)*x1seq
                      + mean(rp2$beta2_mu)*0
                      + mean(rp2$beta3_mu)*0*x1seq
),type = "l", lwd = 2, col = rgb(0,0,0,1), lty = 1)




### - 1 SD
bigN <- nrow(rp2$alpha_mu)
bigN_index <- seq(1,bigN,10)
bigN <- length(bigN_index)
pred1 <- matrix(NA,nrow = bigN,ncol = length(x1seq))
for(i in 1:bigN)  pred1[i,] <- invlogit((rp2$alpha_mu[i]) + rp2$beta1_mu[i]*x1seq + rp2$beta2_mu[i]*(-1) + rp2$beta3_mu[i]*(-1)*x1seq)

low1 <- apply(pred1,2,function(x) quantile(x,0.025))
high1 <- apply(pred1,2,function(x) quantile(x,0.975))
med1 <- apply(pred1,2,function(x) quantile(x,0.5))
mean1 <- apply(pred1,2,function(x) mean(x))

#error_polygon(low1,high1,x1seq[1],x1seq[length(x1seq)],x1seq,
#              col = rgb(0,0.7,0.7,0.3))
points(x1seq,invlogit(mean(rp2$alpha_mu)+mean(rp2$beta1_mu)*x1seq
                      + mean(rp2$beta2_mu)*(-1)
                      + mean(rp2$beta3_mu)*(-1)*x1seq
),type = "l", lwd = 2, col = rgb(0,0.5,0.7,1), lty = 4)





### + 1 SD
bigN <- nrow(rp2$alpha_mu)
bigN_index <- seq(1,bigN,10)
bigN <- length(bigN_index)
pred1 <- matrix(NA,nrow = bigN,ncol = length(x1seq))
for(i in 1:bigN)  pred1[i,] <- invlogit((rp2$alpha_mu[i]) + rp2$beta1_mu[i]*x1seq + rp2$beta2_mu[i]*1 + rp2$beta3_mu[i]*1*x1seq)

low1 <- apply(pred1,2,function(x) quantile(x,0.025))
high1 <- apply(pred1,2,function(x) quantile(x,0.975))
med1 <- apply(pred1,2,function(x) quantile(x,0.5))
mean1 <- apply(pred1,2,function(x) mean(x))

#error_polygon(low1,high1,x1seq[1],x1seq[length(x1seq)],x1seq,
#              col = rgb(0,0.7,0.7,0.3))
points(x1seq,invlogit(mean(rp2$alpha_mu)+mean(rp2$beta1_mu)*x1seq
                      + mean(rp2$beta2_mu)*1
                      + mean(rp2$beta3_mu)*1*x1seq
),type = "l", lwd = 2, col = rgb(0.7,0,0.5,1), lty = 2)


legend("topleft", c(expression("mean"~italic("S")), expression("mean"~italic("S")~"+1 SD"), expression("mean"~italic("S")~"-1 SD")),
       col = c("black", rgb(0.7,0,0.5,1), rgb(0,0.5,0.7,1)), lty = c(1,2,4), lwd = 2, cex = 0.8, box.col = NA, bg = NA)


dev.off()


#########################################################
###
### Fig 3c


#### T range on X

windows(3.5,3.5)
#png("Revision_1/hierarch_reg_mu_T_S_TS_sx+8r.png", units = "in", width = 3.1, height = 3.1, res = 600)
#pdf("hierarch_reg_mu_T_S_TS_sx+8r.pdf", width = 3.1, height = 3.1)

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



sizes <- (dat2$n)^0.11-0.45
par(mar = c(4,4,1,1), las = 1, mgp = c(2.5,0.75,0))
plot(0,0, pch = 21, col = NA, bg = rgb(0,0,0,0.3), cex = sizes, xlab = expression("standardised seasonality"), ylab = "P arag. occurrences",
     xaxs = "i", xlim = c(min(x2)-0.1,max(x2)+0.1), type = "n", ylim = c(0,1))
x1seq <- seq(min(x2),max(x2),0.05)



### Mean
bigN <- nrow(rp2$alpha_mu)
bigN_index <- seq(1,bigN,1)
bigN <- length(bigN_index)
pred1 <- matrix(NA,nrow = bigN,ncol = length(x1seq))
for(i in 1:bigN)  pred1[i,] <- invlogit((rp2$alpha_mu[i]) + rp2$beta1_mu[i]*0 + rp2$beta2_mu[i]*x1seq + rp2$beta3_mu[i]*0*x1seq)

low1 <- apply(pred1,2,function(x) quantile(x,0.025))
high1 <- apply(pred1,2,function(x) quantile(x,0.975))
med1 <- apply(pred1,2,function(x) quantile(x,0.5))
mean1 <- apply(pred1,2,function(x) mean(x))

error_polygon(low1,high1,x1seq[1],x1seq[length(x1seq)],x1seq,
              col = rgb(0,0,0,0.25))
points(x1seq,invlogit(mean(rp2$alpha_mu)+mean(rp2$beta1_mu)*0
                      + mean(rp2$beta2_mu)*x1seq
                      + mean(rp2$beta3_mu)*0*x1seq
),type = "l", lwd = 2, col = rgb(0,0,0,1), lty = 1)




### -1 SD
bigN <- nrow(rp2$alpha_mu)
bigN_index <- seq(1,bigN,10)
bigN <- length(bigN_index)
pred1 <- matrix(NA,nrow = bigN,ncol = length(x1seq))
for(i in 1:bigN)  pred1[i,] <- invlogit((rp2$alpha_mu[i]) + rp2$beta1_mu[i]*(-1) + rp2$beta2_mu[i]*x1seq + 
                                          rp2$beta3_mu[i]*(-1)*x1seq)

low1 <- apply(pred1,2,function(x) quantile(x,0.025))
high1 <- apply(pred1,2,function(x) quantile(x,0.975))
med1 <- apply(pred1,2,function(x) quantile(x,0.5))
mean1 <- apply(pred1,2,function(x) mean(x))

#error_polygon(low1,high1,x1seq[1],x1seq[length(x1seq)],x1seq,
#              col = rgb(0,0.7,0.7,0.3))
points(x1seq,invlogit(mean(rp2$alpha_mu)+mean(rp2$beta1_mu)*(-1)
                      + mean(rp2$beta2_mu)*x1seq
                      + mean(rp2$beta3_mu)*(-1)*x1seq
),type = "l", lwd = 2, col = rgb(0,0.5,0.7,1), lty = 4)





### +1 SD
bigN <- nrow(rp2$alpha_mu)
bigN_index <- seq(1,bigN,10)
bigN <- length(bigN_index)
pred1 <- matrix(NA,nrow = bigN,ncol = length(x1seq))
for(i in 1:bigN)  pred1[i,] <- invlogit((rp2$alpha_mu[i]) + rp2$beta1_mu[i]*1 + rp2$beta2_mu[i]*x1seq + 
                                          rp2$beta3_mu[i]*1*x1seq)

low1 <- apply(pred1,2,function(x) quantile(x,0.025))
high1 <- apply(pred1,2,function(x) quantile(x,0.975))
med1 <- apply(pred1,2,function(x) quantile(x,0.5))
mean1 <- apply(pred1,2,function(x) mean(x))

#error_polygon(low1,high1,x1seq[1],x1seq[length(x1seq)],x1seq,
#              col = rgb(0,0.7,0.7,0.3))
points(x1seq,invlogit(mean(rp2$alpha_mu)+mean(rp2$beta1_mu)*1
                      + mean(rp2$beta2_mu)*x1seq
                      + mean(rp2$beta3_mu)*1*x1seq
),type = "l", lwd = 2, col =  rgb(0.7,0,0.5,1), lty = 2)

legend("topleft", c(expression("mean"~italic("T")), expression("mean"~italic("T")~"+1 SD"), expression("mean"~italic("T")~"-1 SD")),
       col = c("black", rgb(0.7,0,0.5,1), rgb(0,0.5,0.7,1)), lty = c(1,2,4), lwd = 2, cex = 0.8, box.col = NA, bg = NA)

dev.off()