############################################### # PlotFullModelFits.R # Sean Connolly # Cleaned up 27 March 2025 # ----------------------------------------- # This code loads and plots the best-fitting models shown in ED Figures 3-5. ################################################ rm(list = ls()) #### Set up #### # Load required libraries library("glmmTMB") library("latex2exp") library("scales") # Set working directory as appropriate # Load the dataset orig_data = read.csv("Database_v7_061219_fixNA.csv") bleaching_data = as.data.frame(orig_data) mortality_data = as.data.frame(orig_data) # Bleaching/mortality severity code -1 = % unknown # Convert -1 to NA in the dataset bleaching_data$BLEACHING_SEVERITY_CODE_FIXED[which(orig_data$BLEACHING_SEVERITY_CODE_FIXED == -1)] = NA mortality_data$MORTALITY_SEVERITY_CODE_FIXED[which(orig_data$MORTALITY_SEVERITY_CODE_FIXED == -1)] = NA # Keep only the data for which the bleaching/mortality severity is known bleaching_data = bleaching_data[which(is.na(bleaching_data$BLEACHING_SEVERITY_CODE_FIXED)==FALSE),] mortality_data = mortality_data[which(is.na(mortality_data$MORTALITY_SEVERITY_CODE_FIXED)==FALSE),] # Create binomial bleaching/mortality variables (0:unbleached, 1:bleached) # (TRUE if bleaching score >= threshold, FALSE otherwise) bleaching_data$binom.bleach1 = as.numeric(bleaching_data$BLEACHING_SEVERITY_CODE_FIXED >= 2) bleaching_data$binom.bleach2 = as.numeric(bleaching_data$BLEACHING_SEVERITY_CODE_FIXED >= 3) mortality_data$binom.mort1 = as.numeric(mortality_data$MORTALITY_SEVERITY_CODE_FIXED >= 2) #### SET VARIABLES BY COMMENTING OR UNCOMMENTING THE APPROPRIATE LINES #### data.of.interest = bleaching_data # Specify bleaching_data for bleaching model, mortality_data for mortality model which1 = "Bleaching" # which1 is used to paste a y axis label on plots (specify either "Bleaching" or "Mortality") # data.of.interest = mortality_data # Specify bleaching_data for bleaching model, mortality_data for mortality model # which1 = "Mortality" # which1 is used to paste a y axis label on plots (specify either "Bleaching" or "Mortality") load("Bleach1ModelsFinal.Rdata") # load("Bleach2ModelsFinal.Rdata") # load("MortModelsFinal.Rdata") model <- DHW.R.BY.minusR.BY # Bleach1 best model # model <- DHW.R.BY.minusR.BY.minusR.minusBY # Bleach2 best model # model <- DHW.R.BY # Mort 1 best AIC model # model <- DHW.R.BY.minusDHW.R.BY # Best mortality model that gives no convergence # warnings (Note the # top 3 models in model selection give # convergence errors in the current [as of # 27 March 2025] version of glmmTMB running # on R 4.4.2. These convergence errrors were # not reported when the original analyses were # done in 2021.) ylabel = "Probability of Bleaching (>10% corals affected)" # ylabel = "Probability of Severe Bleaching (>50% corals affected)" # ylabel = "Probability of Mortality (>10% corals affected)" ################################################################## #### Plot Model predictions #### # Specify data for predictions pred = expand.grid(DHW_TO_USE = seq(min(data.of.interest$DHW_TO_USE), max(data.of.interest$DHW_TO_USE), len = 500), REGION = factor(c('AP', 'CA', 'IM'), levels = c('AP', 'CA', 'IM')), BLEACHING_YEAR = factor(c('2014-2015', '2015-2016', '2016-2017'), levels = c('2014-2015', '2015-2016', '2016-2017'))) # Create model matrix mm = model.matrix(lme4::nobars(formula(model))[-2], pred) # Calculate model predictions (fixed effects only) model.pred = mm%*%(fixef(model)$cond) # Calculate 95% confidence intervals pvar = diag(mm %*% tcrossprod(vcov(model)$cond, mm)) pred$lowCI = 1/(1+exp(-(model.pred-1.96*sqrt(pvar)))) pred$highCI = 1/(1+exp(-(model.pred+1.96*sqrt(pvar)))) # Plot windows(height = 9) par(mfrow = c(3,1), cex.main = 1.2, omi = c(0.5,0.5,0.2,0.2), mar = c(2,2,4,1)) plot(pred$DHW_TO_USE[which(pred$REGION == "AP" & pred$BLEACHING_YEAR == "2014-2015" & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'AP' & data.of.interest$BLEACHING_YEAR == '2014-2015')]))], 1/(1+exp(-(model.pred[which(pred$REGION == "AP" & pred$BLEACHING_YEAR == "2014-2015" & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'AP' & data.of.interest$BLEACHING_YEAR == '2014-2015')]))]))), type = "l", col = "red", ylab = "", xlab = "", ylim = c(0,1), main = "Asia-Pacific", xlim = c(0,30),cex.axis=2,cex.main=2) polygon(c(rev(pred$DHW_TO_USE[which(pred$REGION == 'AP' & pred$BLEACHING_YEAR == '2014-2015' & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'AP' & data.of.interest$BLEACHING_YEAR == '2014-2015')]))]), pred$DHW_TO_USE[which(pred$REGION == 'AP' & pred$BLEACHING_YEAR == '2014-2015' & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'AP' & data.of.interest$BLEACHING_YEAR == '2014-2015')]))]), c(rev(pred$highCI[which(pred$REGION == 'AP'& pred$BLEACHING_YEAR == '2014-2015' & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'AP' & data.of.interest$BLEACHING_YEAR == '2014-2015')]))]), pred$lowCI[which(pred$REGION == 'AP'& pred$BLEACHING_YEAR == '2014-2015' & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'AP' & data.of.interest$BLEACHING_YEAR == '2014-2015')]))]), col = alpha('red', 0.2), border = NA) polygon(c(rev(pred$DHW_TO_USE[which(pred$REGION == 'AP' & pred$BLEACHING_YEAR == '2015-2016' & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'AP' & data.of.interest$BLEACHING_YEAR == '2015-2016')]))]), pred$DHW_TO_USE[which(pred$REGION == 'AP' & pred$BLEACHING_YEAR == '2015-2016' & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'AP' & data.of.interest$BLEACHING_YEAR == '2015-2016')]))]), c(rev(pred$highCI[which(pred$REGION == 'AP'& pred$BLEACHING_YEAR == '2015-2016' & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'AP' & data.of.interest$BLEACHING_YEAR == '2015-2016')]))]), pred$lowCI[which(pred$REGION == 'AP'& pred$BLEACHING_YEAR == '2015-2016' & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'AP' & data.of.interest$BLEACHING_YEAR == '2015-2016')]))]), col = alpha('blue', 0.2), border = NA) polygon(c(rev(pred$DHW_TO_USE[which(pred$REGION == 'AP' & pred$BLEACHING_YEAR == '2016-2017' & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'AP' & data.of.interest$BLEACHING_YEAR == '2016-2017')]))]), pred$DHW_TO_USE[which(pred$REGION == 'AP' & pred$BLEACHING_YEAR == '2016-2017' & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'AP' & data.of.interest$BLEACHING_YEAR == '2016-2017')]))]), c(rev(pred$highCI[which(pred$REGION == 'AP'& pred$BLEACHING_YEAR == '2016-2017' & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'AP' & data.of.interest$BLEACHING_YEAR == '2016-2017')]))]), pred$lowCI[which(pred$REGION == 'AP'& pred$BLEACHING_YEAR == '2016-2017' & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'AP' & data.of.interest$BLEACHING_YEAR == '2016-2017')]))]), col = alpha('forestgreen', 0.2), border = NA) lines(pred$DHW_TO_USE[which(pred$REGION == "AP" & pred$BLEACHING_YEAR == "2014-2015" & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'AP' & data.of.interest$BLEACHING_YEAR == '2014-2015')]))], 1/(1+exp(-(model.pred[which(pred$REGION == "AP" & pred$BLEACHING_YEAR == "2014-2015" & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'AP' & data.of.interest$BLEACHING_YEAR == '2014-2015')]))]))), col = "red", lwd = 2) lines(pred$DHW_TO_USE[which(pred$REGION == "AP" & pred$BLEACHING_YEAR == "2015-2016" & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'AP' & data.of.interest$BLEACHING_YEAR == '2015-2016')]))], 1/(1+exp(-(model.pred[which(pred$REGION == "AP" & pred$BLEACHING_YEAR == "2015-2016" & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'AP' & data.of.interest$BLEACHING_YEAR == '2015-2016')]))]))), col = "blue", lwd = 2) lines(pred$DHW_TO_USE[which(pred$REGION == "AP" & pred$BLEACHING_YEAR == "2016-2017" & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'AP' & data.of.interest$BLEACHING_YEAR == '2016-2017')]))], 1/(1+exp(-(model.pred[which(pred$REGION == "AP" & pred$BLEACHING_YEAR == "2016-2017" & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'AP' & data.of.interest$BLEACHING_YEAR == '2016-2017')]))]))), col = "forestgreen", lwd = 2) text(0,0.93,labels="a",cex=2) plot(pred$DHW_TO_USE[which(pred$REGION == "CA" & pred$BLEACHING_YEAR == "2014-2015" & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'CA' & data.of.interest$BLEACHING_YEAR == '2014-2015')]))], 1/(1+exp(-(model.pred[which(pred$REGION == "CA" & pred$BLEACHING_YEAR == "2014-2015" & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'CA' & data.of.interest$BLEACHING_YEAR == '2014-2015')]))]))), type = "l", col = "red", ylab = "", xlab = "", ylim = c(0,1), main = "Caribbean/Atlantic", xlim = c(0,30),cex.lab = 2,cex.main=2,cex.axis=2) polygon(c(rev(pred$DHW_TO_USE[which(pred$REGION == 'CA' & pred$BLEACHING_YEAR == '2014-2015' & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'CA' & data.of.interest$BLEACHING_YEAR == '2014-2015')]))]), pred$DHW_TO_USE[which(pred$REGION == 'CA' & pred$BLEACHING_YEAR == '2014-2015' & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'CA' & data.of.interest$BLEACHING_YEAR == '2014-2015')]))]), c(rev(pred$highCI[which(pred$REGION == 'CA'& pred$BLEACHING_YEAR == '2014-2015' & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'CA' & data.of.interest$BLEACHING_YEAR == '2014-2015')]))]), pred$lowCI[which(pred$REGION == 'CA'& pred$BLEACHING_YEAR == '2014-2015' & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'CA' & data.of.interest$BLEACHING_YEAR == '2014-2015')]))]), col = alpha('red', 0.2), border = NA) polygon(c(rev(pred$DHW_TO_USE[which(pred$REGION == 'CA' & pred$BLEACHING_YEAR == '2015-2016' & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'CA' & data.of.interest$BLEACHING_YEAR == '2015-2016')]))]), pred$DHW_TO_USE[which(pred$REGION == 'CA' & pred$BLEACHING_YEAR == '2015-2016' & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'CA' & data.of.interest$BLEACHING_YEAR == '2015-2016')]))]), c(rev(pred$highCI[which(pred$REGION == 'CA'& pred$BLEACHING_YEAR == '2015-2016' & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'CA' & data.of.interest$BLEACHING_YEAR == '2015-2016')]))]), pred$lowCI[which(pred$REGION == 'CA'& pred$BLEACHING_YEAR == '2015-2016' & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'CA' & data.of.interest$BLEACHING_YEAR == '2015-2016')]))]), col = alpha('blue', 0.2), border = NA) polygon(c(rev(pred$DHW_TO_USE[which(pred$REGION == 'CA' & pred$BLEACHING_YEAR == '2016-2017' & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'CA' & data.of.interest$BLEACHING_YEAR == '2016-2017')]))]), pred$DHW_TO_USE[which(pred$REGION == 'CA' & pred$BLEACHING_YEAR == '2016-2017' & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'CA' & data.of.interest$BLEACHING_YEAR == '2016-2017')]))]), c(rev(pred$highCI[which(pred$REGION == 'CA'& pred$BLEACHING_YEAR == '2016-2017' & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'CA' & data.of.interest$BLEACHING_YEAR == '2016-2017')]))]), pred$lowCI[which(pred$REGION == 'CA'& pred$BLEACHING_YEAR == '2016-2017' & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'CA' & data.of.interest$BLEACHING_YEAR == '2016-2017')]))]), col = alpha('forestgreen', 0.2), border = NA) lines(pred$DHW_TO_USE[which(pred$REGION == "CA" & pred$BLEACHING_YEAR == "2014-2015" & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'CA' & data.of.interest$BLEACHING_YEAR == '2014-2015')]))], 1/(1+exp(-(model.pred[which(pred$REGION == "CA" & pred$BLEACHING_YEAR == "2014-2015" & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'CA' & data.of.interest$BLEACHING_YEAR == '2014-2015')]))]))), col = "red", lwd = 2) lines(pred$DHW_TO_USE[which(pred$REGION == "CA" & pred$BLEACHING_YEAR == "2015-2016" & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'CA' & data.of.interest$BLEACHING_YEAR == '2015-2016')]))], 1/(1+exp(-(model.pred[which(pred$REGION == "CA" & pred$BLEACHING_YEAR == "2015-2016" & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'CA' & data.of.interest$BLEACHING_YEAR == '2015-2016')]))]))), col = "blue", lwd = 2) lines(pred$DHW_TO_USE[which(pred$REGION == "CA" & pred$BLEACHING_YEAR == "2016-2017" & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'CA' & data.of.interest$BLEACHING_YEAR == '2016-2017')]))], 1/(1+exp(-(model.pred[which(pred$REGION == "CA" & pred$BLEACHING_YEAR == "2016-2017" & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'CA' & data.of.interest$BLEACHING_YEAR == '2016-2017')]))]))), col = "forestgreen", lwd = 2) legend(20, 0.8, c("2014-2015", "2015-2016", "2016-2017"), col = c("red", "blue", "forestgreen"), pch = 16, bty = "n",cex=2) text(0,0.93,labels="b",cex=2) plot(pred$DHW_TO_USE[which(pred$REGION == "IM" & pred$BLEACHING_YEAR == "2014-2015" & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'IM' & data.of.interest$BLEACHING_YEAR == '2014-2015')]))], 1/(1+exp(-(model.pred[which(pred$REGION == "IM" & pred$BLEACHING_YEAR == "2014-2015" & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'IM' & data.of.interest$BLEACHING_YEAR == '2014-2015')]))]))), type = "l", col = "red", ylab = "", xlab = "", ylim = c(0,1), main = "Indian Ocean/Middle East", xlim = c(0,30),cex.main=2,cex.axis=2,cex.lab=2) polygon(c(rev(pred$DHW_TO_USE[which(pred$REGION == 'IM' & pred$BLEACHING_YEAR == '2014-2015' & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'IM' & data.of.interest$BLEACHING_YEAR == '2014-2015')]))]), pred$DHW_TO_USE[which(pred$REGION == 'IM' & pred$BLEACHING_YEAR == '2014-2015' & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'IM' & data.of.interest$BLEACHING_YEAR == '2014-2015')]))]), c(rev(pred$highCI[which(pred$REGION == 'IM'& pred$BLEACHING_YEAR == '2014-2015' & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'IM' & data.of.interest$BLEACHING_YEAR == '2014-2015')]))]), pred$lowCI[which(pred$REGION == 'IM'& pred$BLEACHING_YEAR == '2014-2015' & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'IM' & data.of.interest$BLEACHING_YEAR == '2014-2015')]))]), col = alpha('red', 0.2), border = NA) polygon(c(rev(pred$DHW_TO_USE[which(pred$REGION == 'IM' & pred$BLEACHING_YEAR == '2015-2016' & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'IM' & data.of.interest$BLEACHING_YEAR == '2015-2016')]))]), pred$DHW_TO_USE[which(pred$REGION == 'IM' & pred$BLEACHING_YEAR == '2015-2016' & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'IM' & data.of.interest$BLEACHING_YEAR == '2015-2016')]))]), c(rev(pred$highCI[which(pred$REGION == 'IM'& pred$BLEACHING_YEAR == '2015-2016' & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'IM' & data.of.interest$BLEACHING_YEAR == '2015-2016')]))]), pred$lowCI[which(pred$REGION == 'IM'& pred$BLEACHING_YEAR == '2015-2016' & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'IM' & data.of.interest$BLEACHING_YEAR == '2015-2016')]))]), col = alpha('blue', 0.2), border = NA) polygon(c(rev(pred$DHW_TO_USE[which(pred$REGION == 'IM' & pred$BLEACHING_YEAR == '2016-2017' & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'IM' & data.of.interest$BLEACHING_YEAR == '2016-2017')]))]), pred$DHW_TO_USE[which(pred$REGION == 'IM' & pred$BLEACHING_YEAR == '2016-2017' & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'IM' & data.of.interest$BLEACHING_YEAR == '2016-2017')]))]), c(rev(pred$highCI[which(pred$REGION == 'IM'& pred$BLEACHING_YEAR == '2016-2017' & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'IM' & data.of.interest$BLEACHING_YEAR == '2016-2017')]))]), pred$lowCI[which(pred$REGION == 'IM'& pred$BLEACHING_YEAR == '2016-2017' & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'IM' & data.of.interest$BLEACHING_YEAR == '2016-2017')]))]), col = alpha('forestgreen', 0.2), border = NA) lines(pred$DHW_TO_USE[which(pred$REGION == "IM" & pred$BLEACHING_YEAR == "2014-2015" & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'IM' & data.of.interest$BLEACHING_YEAR == '2014-2015')]))], 1/(1+exp(-(model.pred[which(pred$REGION == "IM" & pred$BLEACHING_YEAR == "2014-2015" & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'IM' & data.of.interest$BLEACHING_YEAR == '2014-2015')]))]))), col = "red", lwd = 2) lines(pred$DHW_TO_USE[which(pred$REGION == "IM" & pred$BLEACHING_YEAR == "2015-2016" & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'IM' & data.of.interest$BLEACHING_YEAR == '2015-2016')]))], 1/(1+exp(-(model.pred[which(pred$REGION == "IM" & pred$BLEACHING_YEAR == "2015-2016" & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'IM' & data.of.interest$BLEACHING_YEAR == '2015-2016')]))]))), col = "blue", lwd = 2) lines(pred$DHW_TO_USE[which(pred$REGION == "IM" & pred$BLEACHING_YEAR == "2016-2017" & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'IM' & data.of.interest$BLEACHING_YEAR == '2016-2017')]))], 1/(1+exp(-(model.pred[which(pred$REGION == "IM" & pred$BLEACHING_YEAR == "2016-2017" & pred$DHW_TO_USE <= max(data.of.interest$DHW_TO_USE[which(data.of.interest$REGION == 'IM' & data.of.interest$BLEACHING_YEAR == '2016-2017')]))]))), col = "forestgreen", lwd = 2) text(0,0.93,labels="c",cex=2) mtext(TeX("Degree Heating Weeks ($^o$C-weeks)"),side=1,adj=0.5,line=2,cex=1.5,outer=TRUE) mtext(ylabel,side=2,adj=0.5,line=1,cex=1.5,outer=TRUE)