R Code for GCB-19-1650.R1

rm(list = ls())  # Remove all objects from memory

1. Influence of most intense ENSO events on the coral cover in the ETP

library(ggplot2)
library(gridExtra)

ENSO_colours<-c("blue", "red") # Type of anomaly
my_colours<-c("#feb24c93", "#fb6a4a93", "#9e9ac893", "#a6d96a93") # Thermal regime

theme_set (theme_classic() + theme(panel.grid.major = element_blank(),
                              panel.grid.minor = element_blank(), 
                              axis.line = element_line(colour = "black"),
                              legend.position="none",
                              axis.text.x = element_text(angle = 90, vjust = 0.5),
                              plot.title = element_text(size=12, face="bold"),
                              #panel.border = element_rect(colour = "black", fill=NA, size=1)
                              panel.border = element_blank()
                              ))

ENSO intensity in region 3

Figure 2a

# Plot ENOS 1971 - 2014
  ENOS<-read.table(file="http://www.cpc.ncep.noaa.gov/data/indices/ersst5.nino.mth.81-10.ascii",header=TRUE)
  ENOS1982_2014 <-subset(ENOS, YR>1968 & YR<2015)   # subset from 1970-2014
  ENOS1982_2014$Date<-paste(ENOS1982_2014$YR, ENOS1982_2014$MON, sep = "-" )
  ENOS1982_2014$Date<-paste(ENOS1982_2014$Date, "15", sep = "-" )
  ENOS1982_2014$Date<-as.Date(ENOS1982_2014$Date, format = "%Y-%m-%d")
  ENOS1982_2014$Anomaly <- factor(ifelse(ENOS1982_2014$ANOM.3 >=0, "Positive", "Negative"))
  
  interp <- approx(ENOS1982_2014$Date, ENOS1982_2014$ANOM.3, n=10000)
  DataENSO <- data.frame(Date=interp$x, Anomaly3=interp$y)
  DataENSO$Anomaly <- factor(ifelse(DataENSO$Anomaly3>=0, "Positive", "Negative"))
      
ENSO<-ggplot(data=DataENSO, aes(x=Date, y=Anomaly3))+
  scale_x_continuous("",expand=c(0, 0))+
  scale_y_continuous("SST anomaly (ºC) in Niño.3",expand=c(0, 0))+
  geom_area(aes(fill=Anomaly), alpha=0.7)+
  scale_fill_manual(values=ENSO_colours) + labs(title="a")+
  geom_hline(yintercept=0, linetype="solid", color="black")+
  geom_line(position = "identity", size=0.3)+
  theme(axis.text.x=element_blank(),
        axis.title.x=element_blank(),
        axis.ticks.x = element_blank())+
  annotate("text", x = 1900, y = 2.0, label = "El Niño \n72-73", size=3)+
  annotate("text", x = 5600, y = 2.0, label = "El Niño \n82-83", size=3)+
  annotate("text", x = 11100, y = 2.0, label = "El Niño \n97-98", size=3)
ENSO

# ENSOb<-ggplot(data=ENOS1982_2014, aes(x=Date, y=ANOM.3))+
#   geom_hline(yintercept=0, linetype="solid", color="black")+
#   scale_x_date("", limit=c(as.Date("1969-01-15"), as.Date("2014-12-31")),
#                breaks= "12 months",
#                      expand=c(0, 0))+
#   scale_y_continuous("SST anomaly (ºC) in Niño.3",expand=c(0, 0))+
#   geom_area(aes(fill=Anomaly))+
#   geom_line(position = "identity")+
#   scale_fill_manual(values=ENSO_colours) + labs(title="(a)")
# ENSOb

2. DHW and coral cover change

# Load packages 
library(tidyverse)
library(dplyr)
library(data.table)

2.1 Yearly maximum DHW experienced by the best-represented coral reefs sites from 1982 to 2014

  • Select and transform DHW data
# Select DHW data file: maxDHW.csv
  dhw <- read.csv("Data/maxDHW.csv", header=TRUE)

# Exclude years 1981, and 2017-2019
  dhw <- dhw[!dhw$Year == 1981, ]

# Create groups for plot
  ts.Cano_island <- subset(dhw, dhw$Site == "Cano_island")     
  ts.Uva_island <- subset(dhw, dhw$Site == "Uva_Island")          
  ts.Gorgona_island <- subset(dhw, dhw$Site == "Gorgona_Island")          
  ts.Darwin_island <- subset(dhw, dhw$Site == "Darwin_North_Anchorage")          
  ts.Bahia_Banderas <- subset(dhw, dhw$Site == "Bahia_Banderas")
  
  dhw.best.represented<-rbind(ts.Cano_island, ts.Uva_island, ts.Gorgona_island, 
                              ts.Darwin_island, ts.Bahia_Banderas)

FIGURE 3a

  • Plot max DHW in each location
dhw.best.represented$Thermal_regime<-factor(dhw.best.represented$Site, levels =
                              c("Cano_island", "Uva_Island", "Gorgona_Island", 
                                "Darwin_North_Anchorage", 
                                "Bahia_Banderas"))


# Figure 3a
Figure3a<-ggplot(dhw.best.represented, aes(x=Year, y=maxDHW,
                                 colour=Thermal_regime)) + theme_bw() +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(), 
        axis.text.x = element_text(angle = 90, vjust = 0.5),
        legend.background = element_blank(),
        legend.position = c(0.77,0.77),
        legend.title = element_blank(),
        legend.text = element_text(size=8))+
  ggtitle("a") +
  scale_x_continuous(limits = c(1982,2014),
                     breaks = seq(1982,2014, 1),
                     expand = c(0.05 , 0.05),
                      "Years")+
  scale_y_continuous(limits = c(0, 28),
                     breaks = seq(0,28, 4),
                    "Degree heating weeks (°C-weeks)")+
  
  annotate("rect", xmin = 1982, xmax = 1983.5, ymin = 0, ymax = 28, fill="#fee8c8", alpha = .3)+
  annotate("rect", xmin = 1997, xmax = 1998.5, ymin = 0, ymax = 28, fill="#fee8c8", alpha = .3)+
  geom_hline(yintercept=4, linetype=2)+
  geom_hline(yintercept = 8, linetype=2)+
  geom_line() +
  scale_colour_manual(values=c("#8c2d04","#feb24c", "#fb6a4a", "#9e9ac8", "#a6d96a"),
                    labels=c( "Thermally stable - Caño Island", 
                              "Thermally stable - Uva Island", 
                              "Tropical upwelling - Gorgona Island",
                              "Equatorial upwelling - Darwin Island",
                              "Seasonal - Bahía Banderas"))

2.2 Relationship of coral cover rate of change and thermal stress

  • The DHW index measures the accumulated thermal stress above a bleaching threshold for 12 consecutive weeks.
  • At each site, the maximum monthly mean sea surface # temperature (MMM) is calculated from a long-term satellite climatology dataset.
  • Any temperature surpassing the MMM by 1ºC or more is added for 12 consecutive weeks, because temperatures ~1ºC above the MMM are known to cause coral bleaching (Liu et al., 2006).
  • For example, if the water temperature during four consecutive weeks exceeds the MMM in 1.3°C, # 1.1°C, 1.2°C, and 1.4°C, then, the site accumulates 5°C-weeks when these anomalies are summed.

We tested the hypothesis that the strength of the Log_annual_rate_of_change (response variable) is a function of the maxDHW experienced. We had measured the rate of change from four Thermal regimes, and fitted Thermal regimes as a random intercept, and estimate a common slope (change in the rate of change) for maxDHW across all Thermal regimes by fitting it as a fixed effect.

## Load data
  roc.cover <- read.csv("Data/RoC_DHW.csv", header=TRUE)  ### Select file:  RoC_DWH.csv

# Load packages
  library(lme4)
  library(lmerTest)
  library(MuMIn)

Model 6: coral cover change and DHW

model6 <- lmerTest::lmer(Log_annual_rate_of_change ~ maxDHW + (1|Thermal_regime), data=roc.cover)
step(model6)
## Backward reduced random-effect table:
## 
##                      Eliminated npar  logLik    AIC    LRT Df Pr(>Chisq)
## <none>                             4 -17.738 43.476                     
## (1 | Thermal_regime)          0    3 -27.772 61.545 20.069  1  7.469e-06
##                         
## <none>                  
## (1 | Thermal_regime) ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##        Eliminated  Sum Sq Mean Sq NumDF  DenDF F value   Pr(>F)   
## maxDHW          0 0.65054 0.65054     1 128.59   9.873 0.002082 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## Log_annual_rate_of_change ~ maxDHW + (1 | Thermal_regime)
summary(model6) 
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Log_annual_rate_of_change ~ maxDHW + (1 | Thermal_regime)
##    Data: roc.cover
## 
## REML criterion at convergence: 35.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.7962 -0.1472  0.1167  0.4152  1.8621 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  Thermal_regime (Intercept) 0.02501  0.1581  
##  Residual                   0.06589  0.2567  
## Number of obs: 133, groups:  Thermal_regime, 4
## 
## Fixed effects:
##               Estimate Std. Error         df t value Pr(>|t|)   
## (Intercept)  -0.085763   0.084256   3.155381  -1.018  0.38032   
## maxDHW       -0.016316   0.005193 128.594887  -3.142  0.00208 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##        (Intr)
## maxDHW -0.165
r.squaredGLMM(model6)   #  returns the marginal and the conditional R2
##             R2m       R2c
## [1,] 0.05249512 0.3131505
# marginal R2  describes the proportion of variance explained by the fixed factor(s) alone:
# conditional R2 describes the proportion of variance explained by both the fixed and random factors
# Nakagawa, S. & Schielzeth, H. (2013).

FIGURE 3b: Relationship of coral cover rate of change and thermal stress

roc.cover$ENSO_Year<-factor(roc.cover$ENSO_Year, levels =
                              c("No", "1982-1984", "1997-1998"))
roc.cover$Thermal_regime<-factor(roc.cover$Thermal_regime, levels =
                              c("Thermally stable", "Equatorial upwelling",
                                "Tropical upwelling", "Seasonal"))

# Figure 3b
Figure3b<-ggplot(roc.cover, aes(x=maxDHW, y=(Log_annual_rate_of_change),
                                 fill=Thermal_regime, shape=ENSO_Year)) + 
  ggtitle("b") +
  scale_x_continuous(limits = c(0,30),
                     breaks = seq(0,30, 5),
                     expand = c(0.05 , 0.05),
                      "Degree heating weeks (°C-weeks)")+
  scale_y_continuous(limits = c(-1.7, 1),
                     breaks = seq(-1.5,1, 0.5),
    "Annual rate of change in coral cover")+
  geom_point(size=3, alpha=0.5) +
  geom_abline(slope = 0, intercept = 0, linetype=2)+
  geom_vline(xintercept=4, linetype=2)+
  geom_vline(xintercept = 8, linetype=2)+
  scale_fill_manual(values=c("#feb24c93", "#9e9ac893", "#fb6a4a93", "#a6d96a93"),
                    labels=c( "Thermally stable", "Equatorial upwelling", "Tropical upwelling", "Seasonal"))+
  scale_shape_manual(values=c(21, 22, 24),
                     labels=c("No ENSO", "1982-1984", "1997-1998"))+
theme_bw() + theme(panel.grid.major = element_blank(),
                  panel.grid.minor = element_blank(),
                  legend.text = element_text(size=8),
                  legend.title = element_blank(),
                  legend.position = c(0.68, 0.6),
                  #legend.position = "right",
                  legend.background = element_blank())
# Figure3b
grid.arrange(Figure3a, Figure3b, nrow = 1, widths = c(3/5, 2/5))

# ggsave(file="Figure_3.svg", dpi = 300, 
# plot=grid.arrange(Figure3a, Figure3b, nrow = 1, widths = c(3/5, 2/5)), width = 12, height = 4, units = "in")

FIGURE 3 Heat stress and its effect on coral cover across the ETP. a The sequence of yearly maximum DHW experienced by the best-represented individual coral reefs sites from 1982 to 2014. When cumulative thermal stress reaches 4°C-weeks and 8°C-weeks or higher (horizontal dotted lines), bleaching is likely, as well as coral mortality from thermal stress (Liu et al., 2006), respectively. Although the thresholds were developed in seasonal systems, and are not universally applicable for equatorial regions, they are still useful for separating responses into groups. b The association between the maximum DHW and the coral cover change at each of four thermal regimes (dot colors as Fig. 2). Each point (n = 133) represents the DHW and the coral cover annual rate of change per individual coral reef site. The horizontal dotted line divides the positive (> 0) and the negative (< 0) rate of change caused by coral growth or mortality, respectively. The vertical dotted line indicates the DHW thresholds of 4°C-weeks and 8°C-weeks, respectively. The vertical and horizontal dotted lines create four response quadrants: I) heat stress and positive coral cover change (resilient response), II) no heat stress and positive coral cover change, III) no heat stress and negative coral cover change, and IV) heat stress and negative coral cover change.

Descriptive statistics heat stress

No stress

no.stress <- subset(roc.cover, maxDHW < 4)                        # subset of DHW < 4
length(no.stress$maxDHW)/length(roc.cover$maxDHW)                 # ratio no stress: total              
## [1] 0.8496241
1- length(no.stress$maxDHW)/length(roc.cover$maxDHW)              # ratio stress: total              
## [1] 0.1503759
positive.rate <- subset(no.stress, Log_annual_rate_of_change > 0) # subset of Rate > 0
length(positive.rate$Log_annual_rate_of_change)                     # Number of observations in 
## [1] 48
negative.rate <- subset(no.stress, Log_annual_rate_of_change < 0)   #subset of Rate > 0
length(negative.rate$Log_annual_rate_of_change)                     # Number of observations in 
## [1] 58

No stress thermally stable

no.stress.Thermally_Stable <- subset(no.stress, Thermal_regime == "Thermally stable")   
positive.rate <- subset(no.stress.Thermally_Stable, Log_annual_rate_of_change > 0)  
negative.rate <- subset(no.stress.Thermally_Stable, Log_annual_rate_of_change < 0)  
length(positive.rate$Log_annual_rate_of_change)/length(negative.rate$Log_annual_rate_of_change)  
## [1] 1.2

No stress tropical upwelling

no.stress.upwelling <- subset(no.stress, Thermal_regime == "Tropical upwelling")
positive.rate <- subset(no.stress.upwelling, Log_annual_rate_of_change > 0)  
negative.rate <- subset(no.stress.upwelling, Log_annual_rate_of_change < 0)  
length(positive.rate$Log_annual_rate_of_change)/length(negative.rate$Log_annual_rate_of_change)  
## [1] 1

No stress equatorial upwelling

no.stress.Galapagos_Islands <- subset(no.stress, Thermal_regime == "Equatorial upwelling") 
positive.rate <- subset(no.stress.Galapagos_Islands, Log_annual_rate_of_change > 0)  
negative.rate <- subset(no.stress.Galapagos_Islands, Log_annual_rate_of_change < 0)  
length(positive.rate$Log_annual_rate_of_change)/length(negative.rate$Log_annual_rate_of_change)  
## [1] 1

No stress seasonal

no.stress.seasonal <- subset(no.stress, Thermal_regime == "Seasonal")   
positive.rate <- subset(no.stress.seasonal, Log_annual_rate_of_change > 0)  
negative.rate <- subset(no.stress.seasonal, Log_annual_rate_of_change < 0)  
length(positive.rate$Log_annual_rate_of_change)/length(negative.rate$Log_annual_rate_of_change)   
## [1] 0.07142857

3. Supplementary figures

library(pals)
library(dplyr)
library(reshape2)
library(RColorBrewer)

Data sources

FIGURE S1: Number of regions and countries surveyed per year

# Barplot number of surveys per year
  obs_num = as.data.frame(table(cover$Year,cover$Region))    
  colnames(obs_num) = c("Year", "Region", "Freq")
  wide_obs_num = dcast(obs_num, Region~Year, value = "Freq")  
  df = data.matrix(wide_obs_num[,2:45])   
  as.table(df)     
##   1970 1971 1972 1974 1975 1976 1977 1978 1979 1980 1981 1982 1983 1984
## A    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## B    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## C    0    0    1    0    0    0    0    0    1    1    0    0    1    0
## D    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## E    0    0    0    0    0    0    0    0    0    1    0    1    0    4
## F    0    0    0    2    3    2    0    0    0    0    0    1    0    1
## G    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## H    3   10    2   12    3    5    3    5    2    5    2    3    5    3
## I    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## J    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## K    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## L    0    0    0    0    0    0    0    2    0    0    0    0    0    0
## M    0    0    0    0    0    0    0    0    0    0    0    1    0    0
## N    2    6    0    2    1    0    0    0    0    0    0    1    4    2
## O    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## P    0    0    0    0    0    0    0    0    0    0    0    0    0    0
##   1985 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998
## A    0    0    0    0    0    0    0    0    0    1    0    0    0    0
## B    0    0    5    0    0    0    0    1    0    0    0    0    0    0
## C    0    0    1    2    4    0    0    1    0    0    2    2    0    4
## D    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## E   12    1    2    0    1    0    0   14    0    5    3    6    0    1
## F    0    1    1    0    1    0    0    0    2    0    0    0    0    0
## G    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## H    4    5    2    4    1    1    0    1    1    1    1    0    1    1
## I    0    0    0    0    0    0    2    0    0    0    0    0   16   10
## J    0    0    1    1    1    0    3    1    1    0    0    0    2    2
## K    0    0    0    0    0    0    0    0    0    0    0    0    2    1
## L    0    0    0    0    0    0    0    0    0    1    1    4    4    2
## M    0    0    0    0    0    0    0    1    0    0    0    0    0    0
## N    0    1    5    1    2    0    1    1    1    1    1    1    1    1
## O    0    0    0    0    0    0    0    0    1    1    0    0    0    0
## P    0    0    0    0    0    0    0    0    0    0    0    0    0    0
##   1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012
## A    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## B    0    0    0    1    0    0    0    1    0    0    4    0    0    0
## C    7    4    3    4    6    4    3    6    3    2    3    2    2    3
## D    0    0    0    0    0    0    0    0    0    0    0    4    0    1
## E    1    0    5    0    4    2    2    0    1    0    2    0    0    0
## F    0    0    1    1    0    1    3    6    2    2    0    0    1    0
## G    0    0    0    0    0    0    0    0    0    0    1    0    0    0
## H    0    1    1   14   12    1    1    1    1    0    1    2    0    0
## I    0    0    0    6    3    0    1    1    4    0   17    6    2    0
## J    0    0    0    1    2    0    0    0    0    5    7    0    0    1
## K    1    5    5    3    2    0    0    0    0    0    2    0    0    0
## L    0    0    2    0    1    0    1    3    2    0   17    3    1    1
## M    0    0    0    0    0    1    2   14    8    0    0    0    0    1
## N    1    1    1    1    1    1    0    0    0    0    1    0    0    0
## O    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## P    0    0    0    0    0    2    0    0    0    0    0    0    0    0
##   2013 2014
## A    0    0
## B    1    1
## C    2    1
## D    4    2
## E    0    0
## F    0    0
## G    0    0
## H    0    0
## I    2    2
## J    0    0
## K    0    0
## L    1    3
## M    0    0
## N    0    0
## O    0    0
## P    0    0
  df[ df>1 ] <- 1 
  levels(cover$Region) # How many regions? = how many colours
##  [1] "Clipperton"                 "Cocos_Islands"             
##  [3] "Colombia"                   "Costa_Rica_Central"        
##  [5] "Costa_Rica_Southern"        "Eastern_Galapagos_Islands" 
##  [7] "El_Salvador"                "Gulf_of_Chiriqui"          
##  [9] "Mexican_Central"            "Mexican_Northern"          
## [11] "Mexican_Southern"           "Nicaragua_Papagayo_Zone"   
## [13] "Northern_Galapagos_Islands" "Panama_Bight"              
## [15] "Revillagigedo_Islands"      "Western_Galapagos_Islands"
  Regions <- c(levels(obs_num$Region))
  color.regions = warmcool(17)
  
  # Figure 2a
  par(mfcol=c(2,1), cex.main = 1.5, mar = c(4,4,1,1), cex.lab = 1.5, cex.axis = 1,  las = 1)
  barplot(df,
          ylim=c(0,15),
          col = color.regions,   
          ylab="Number of regions surveyed per year",  
          xlab="Year",
          #lwd = 1,          
          las = 2,           
          axis.lty = 1, 
          cex.names = 0.4,   
          cex.axis = 0.8)
# Legend as a box
  legend( #"bottomright",    
     x = 2, y = 16, 
     xpd = NA,  # or TRUE 
     legend = Regions,          
     fill = color.regions[1:17],   
     inset=c(-0.5,0),
     cex = 0.35,
     text.font=0.5,            
     x.intersp = -0.5,
     y.intersp = 1,
     lwd = 0.1,
     bty="n")

# Figure 2b
obs_num = as.data.frame(table(cover$Year,cover$Country))    # Table per country
colnames(obs_num) = c("Year", "Country", "Freq")
wide_obs_num = dcast(obs_num, Country~Year, value = "Freq")   
dfc = data.matrix(wide_obs_num[,2:45])   
as.table(dfc)       
##   1970 1971 1972 1974 1975 1976 1977 1978 1979 1980 1981 1982 1983 1984
## A    0    0    1    0    0    0    0    0    1    1    0    0    1    0
## B    0    0    0    0    0    0    0    2    0    1    0    1    0    4
## C    0    0    0    2    3    2    0    0    0    0    0    2    0    1
## D    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## E    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## F    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## G    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## H    5   16    2   14    4    5    3    5    2    5    2    4    9    5
##   1985 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998
## A    0    0    1    2    4    0    0    1    0    0    2    2    0    4
## B   12    1    7    0    1    0    0   15    0    6    4   10    4    3
## C    0    1    1    0    1    0    0    1    2    0    0    0    0    0
## D    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## E    0    0    0    0    0    0    0    0    0    1    0    0    0    0
## F    0    0    1    1    1    0    5    1    2    1    0    0   20   13
## G    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## H    4    6    7    5    3    1    1    2    2    2    2    1    2    2
##   1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012
## A    7    4    3    4    6    4    3    6    3    2    3    2    2    3
## B    1    0    7    1    5    2    3    4    3    0   12    7    1    2
## C    0    0    1    1    0    4    5   20   10    2    0    0    1    1
## D    0    0    0    0    0    0    0    0    0    0    1    0    0    0
## E    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## F    1    5    5   10    7    0    1    1    4    5   26    6    2    1
## G    0    0    0    0    0    0    0    0    0    0   11    0    0    0
## H    1    2    2   15   13    2    1    1    1    0    2    2    0    0
##   2013 2014
## A    2    1
## B    6    6
## C    0    0
## D    0    0
## E    0    0
## F    2    2
## G    0    0
## H    0    0
dfc[ dfc>1 ] <- 1 # replace all values greater than one with one
color.countries <- brewer.pal(9, "Oranges")
countries = c(levels(obs_num$Country))
barplot(dfc,
        col = color.countries,   
        ylim=c(0,10),
        ylab="Number of countries surveyed per year",  
        xlab="Year",
        las = 2,           
        axis.lty = 1,      
        cex.names = 0.4,   
        cex.axis = 0.8)
# Legend as a box
legend( #"bottomright",    
   x = 2, y = 10, 
   xpd = NA,  # or TRUE 
   legend = countries,         
   fill = color.countries[1:9],   
   inset=c(-0.5,0),
   cex = 0.35,
   text.font=1,            
   x.intersp = -0.5,
   y.intersp = 1,
   lwd = 0.1,
   bty="n")

FIGURE S1 Regions and countries surveyed per year. (a) Number of regions surveyed per year in the ETP. Since 1970, the number of surveys in all regions increased by year up to a maximum of 10 regions surveyed in 2009. (b) Number of countries surveyed per year (Colombia, Costa Rica, Ecuador, El Salvador, France, Mexico, Nicaragua and Panama). Since 1970, the number of surveys in all the countries reached a maximum of 6 countries in 2009.

Divide data to test robustness of the patterns

FIGURE S2: Variation in the live coral cover for 1970–2014 per thermal regime

Temporal variation in the live coral cover for 1970–2014.

    1. Thermally stable regime.
    1. Tropical upwelling regime.
    1. Equatorial upwelling regime.
    1. Seasonal regime.

The smooth line represents the best fit to the time series with a 95% confidence level interval and a span = 0.3.

library(ggplot2)
library(lemon)

aggr.region_S2<-aggr.region
summary(aggr.region_S2$Thermal_regime)
##     Thermally stable   Tropical upwelling Equatorial upwelling 
##                   75                   75                   25 
##             Seasonal 
##                   28
aggr.region_S2$Thermal_regime<-factor(aggr.region_S2$Thermal_regime, 
                                  levels=c("Thermally stable", "Tropical upwelling",
                                            "Equatorial upwelling", "Seasonal"),
                                  labels = c ("a. Thermally stable", "b. Tropical upwelling",
                                              "c. Equatorial upwelling", "d. Seasonal"))
                                   

FigureS2<-ggplot(aggr.region_S2, aes(x=Year, y=Coral_cover)) + 
  theme_bw() + theme(panel.border = element_blank(),   # set ggplot to white background
                              panel.grid.major = element_blank(),
                              panel.grid.minor = element_blank(), 
                              axis.line = element_line(),
                              strip.background = element_blank(),
                              strip.text = element_text(face="bold",hjust = 0, vjust = -1))+
  geom_point()+ geom_smooth (span = 0.3) +
  scale_x_continuous(limits=c(1970, 2014), 
                     breaks = seq(1970, 2014, by=5),
                     expand = c(0.01, 0.01), "Year") +
  scale_y_continuous(limits=c(-19, 91), expand = c(0, 0), "Live coral cover (%)") +
  facet_wrap(Thermal_regime~., ncol=1, scales = "free_x", strip.position ="top")

FigureS2

#ggsave(file="FigureS2.svg", plot=FigureS2, width=6.69, height=6.69)

FIGURE S2 Temporal variation in the live coral cover for 1970–2014. (a) Thermally stable regime (n = 75). (b) Tropical upwelling regime (n = 75). (c) Equatorial upwelling regime (n = 25). (d) Seasonal regime (n = 28). The smooth line represents the best fit to the time series with a 95% confidence level interval and a span = 0.3.

FigureS2b<-ggplot(aggr.region_S2, aes(x=Year, y=Coral_cover)) + 
  theme_bw() + theme(panel.border = element_blank(),   # set ggplot to white background
                              panel.grid.major = element_blank(),
                              panel.grid.minor = element_blank(), 
                              axis.line = element_line(),
                              strip.background = element_blank(),
                              strip.text = element_text(face="bold",hjust = 0, vjust = -1))+
  geom_smooth (span = 0.3, colour="gray") +  geom_point(aes(colour=Region)) +
  scale_x_continuous(limits=c(1970, 2014), 
                     breaks = seq(1970, 2014, by=5),
                     expand = c(0.01, 0.01), "Year") +
  scale_y_continuous(limits=c(-19, 91), expand = c(0, 0), "Live coral cover (%)") +
  facet_wrap(Thermal_regime~., ncol=1, scales = "free_x", strip.position ="top") +
  theme(legend.position = "right", legend.direction = "vertical" )

FigureS2b

#ggsave(file="FigureS2.svg", plot=FigureS2b, width=6.69, height=6.69)

FIGURE S2 Temporal variation in the live coral cover for 1970–2014. (a) Thermally stable regime (n = 75). (b) Tropical upwelling regime (n = 75). (c) Equatorial upwelling regime (n = 25). (d) Seasonal regime (n = 28). The smooth line represents the best fit to the time series with a 95% confidence level interval and a span = 0.3. Colors represent the regions sampled in each thermal regime.

FIGURE S3: Data only for four best monitoring sites

We found in published studies four comprehensive monitoring datasets for Uva Island, Gorgona Island, Costa Rica, and the Eastern and Northern Galapagos Islands, which correspond to well recognized long-term monitoring programs

# Create groups
Uva_Island <- subset(aggr.site, Site == "Uva_Island")
Gorgona <-subset(aggr.site, Site == "La_Azufrada")
Cano_island <-subset(aggr.site, Site == "Cano_island")
EGI <-subset(aggr.site, Region == "Eastern_Galapagos_Islands")

Monitoring_sites <- rbind(Uva_Island,Gorgona,Cano_island,EGI)

Monitoring_sites$Thermal_regime<-factor(Monitoring_sites$Thermal_regime,
                  levels=c("Thermally stable", "Tropical upwelling",
                           "Equatorial upwelling", "Seasonal"))

Model 7: Only monitored sites

model7 <- lme(Coral_cover ~ -1 + Thermal_regime * Year, random = ~1|Region, data=Monitoring_sites) 
summary(model7)
## Linear mixed-effects model fit by REML
##  Data: Monitoring_sites 
##        AIC      BIC    logLik
##   736.7157 756.7142 -360.3579
## 
## Random effects:
##  Formula: ~1 | Region
##          (Intercept) Residual
## StdDev: 0.0007454146 10.93506
## 
## Fixed effects: Coral_cover ~ -1 + Thermal_regime * Year 
##                                             Value Std.Error DF   t-value
## Thermal_regimeThermally stable          -575.3388  307.9802  1 -1.868103
## Thermal_regimeTropical upwelling         262.8581  520.3591  1  0.505148
## Thermal_regimeEquatorial upwelling      1078.2137  310.2684  1  3.475100
## Year                                       0.2977    0.1547 90  1.923978
## Thermal_regimeTropical upwelling:Year     -0.4036    0.3025 90 -1.334386
## Thermal_regimeEquatorial upwelling:Year   -0.8301    0.2195 90 -3.782741
##                                         p-value
## Thermal_regimeThermally stable           0.3129
## Thermal_regimeTropical upwelling         0.7022
## Thermal_regimeEquatorial upwelling       0.1784
## Year                                     0.0575
## Thermal_regimeTropical upwelling:Year    0.1854
## Thermal_regimeEquatorial upwelling:Year  0.0003
##  Correlation: 
##                                         Thr_Ts Thr_Tu Thr_Eu Year   T_Tu:Y
## Thermal_regimeTropical upwelling         0.000                            
## Thermal_regimeEquatorial upwelling       0.000  0.000                     
## Year                                    -1.000  0.000  0.000              
## Thermal_regimeTropical upwelling:Year    0.511 -0.859  0.000 -0.511       
## Thermal_regimeEquatorial upwelling:Year  0.705  0.000 -0.709 -0.705  0.361
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -1.8628050 -0.7098789 -0.1002282  0.6820665  2.4590209 
## 
## Number of Observations: 96
## Number of Groups: 4
anova(model7)
plot(ranef(model7))    
plot(model7)           
resnorm7 <- resid(model7)
hist(resnorm7, xlab = "Residuals", main = "") 
coef.m7 <- as.data.frame(coef(summary(model7))) 
plot(model7, resid(., type = "p") ~ fitted(.) | Region, abline = 0)
plot(model7, Coral_cover ~ fitted(.) | Region, abline = c(0,1))
plot(model7, Coral_cover ~ fitted(.) | Thermal_regime, abline = c(0,1))

Figure S3: Data only for four best monitoring sites. a) pooled; b) by Location

summary(Monitoring_sites$Region)
##                 Clipperton              Cocos_Islands 
##                          0                          0 
##                   Colombia         Costa_Rica_Central 
##                         22                          0 
##        Costa_Rica_Southern  Eastern_Galapagos_Islands 
##                         13                         28 
##                El_Salvador           Gulf_of_Chiriqui 
##                          0                         33 
##            Mexican_Central           Mexican_Northern 
##                          0                          0 
##           Mexican_Southern    Nicaragua_Papagayo_Zone 
##                          0                          0 
## Northern_Galapagos_Islands               Panama_Bight 
##                          0                          0 
##      Revillagigedo_Islands  Western_Galapagos_Islands 
##                          0                          0
Monitoring_sites$Region<-factor(Monitoring_sites$Region,
                  levels=c("Costa_Rica_Southern", "Gulf_of_Chiriqui",
                           "Colombia", "Eastern_Galapagos_Islands"))
Region.labels <- as_labeller(c(`Costa_Rica_Southern` = "Thermally stable",
                               `Gulf_of_Chiriqui` = "Thermally stable",
                               `Colombia` = "Tropical upwelling", 
                               `Eastern_Galapagos_Islands` = "Equatorial upwelling"))


FigureS3 <-ggplot(Monitoring_sites, aes(x=Year, y=Coral_cover)) +
  geom_boxplot(aes(group=Year), outlier.shape = NA) +
  stat_boxplot(aes(group=Year), geom = 'errorbar')+
  geom_smooth(method = lm, se=FALSE, linetype = "dashed", colour="red")+
  scale_y_continuous("Live coral cover (%)", limits = c(-10, 95), 
                     breaks = seq(0, 90, by=20), expand = c(0,0))+
  scale_x_continuous("", limits = c(1969, 2015),
                     breaks = seq(1970, 2014, by=2), expand = c(0,0))+
  annotate("rect", xmin = 1982, xmax = 1983, ymin = 0, ymax = 80,  alpha = .2, fill="red")+
  annotate("rect", xmin = 1997, xmax = 1998, ymin = 0, ymax = 80,  alpha = .2, fill="red")+
  annotate("rect", xmin = 1972, xmax = 1973, ymin = 0, ymax = 80,  alpha = .2, fill="red")

#FigureS3

FigureS3a<- FigureS3 + geom_point(aes(colour=Thermal_regime), alpha=0.8) +
  theme(legend.position = "bottom") + ggtitle("a") +
    scale_colour_manual(values=my_colours) +
  geom_smooth(span = 0.3, se=T, colour="darkgray")
#FigureS3a

FigureS3b<-FigureS3 + geom_point(aes(colour=Location), alpha=0.8) +
  facet_wrap(~Region, labeller=Region.labels)+
  geom_smooth(span = 0.3, se=F, colour="darkgray") +
  ggtitle("b") + theme(legend.position = "bottom") 
#FigureS3b

FigureS3_all<-grid.arrange(arrangeGrob(FigureS3a,FigureS3b,
                                       heights=c(1.2/3, 1.8/3)))

#ggsave(file="FigS3.svg", plot=FigureS3_all, width=6.69, height=7.69)

FIGURE S3 Temporal variation in the live coral cover for 1970–2014 for the four best-represented monitoring locations with its fitted linear trend (red dashed band). (a) Average live-coral cover trend for the four best-represented monitoring locations. (b) Comprehensive monitoring datasets were found in published studies for Costa Rica Southern (Caño Island), Gulf of Chiriquí (Uva Island), Colombia (Gorgona Island), and Eastern Galapagos Islands (Floreana, La Espanola, Marchena, San Cristobal, Santa Cruz, and Santa Fe).

4. Extra figures / models

Data excluding one thermal regime from model 1

Model 1: All years (1970_2014) removing each thermal regime

  • Remove thermally stable
 model1_Model1_NoThermallyStable <- lme(
    Coral_cover ~ -1 + Thermal_regime * Year, random = ~1|Region,
    data=subset(aggr.region, Thermal_regime!="Thermally stable"))
  summary(model1_Model1_NoThermallyStable)
## Linear mixed-effects model fit by REML
##  Data: subset(aggr.region, Thermal_regime != "Thermally stable") 
##        AIC      BIC    logLik
##   1121.751 1144.183 -552.8754
## 
## Random effects:
##  Formula: ~1 | Region
##         (Intercept) Residual
## StdDev:    14.39849 18.48283
## 
## Fixed effects: Coral_cover ~ -1 + Thermal_regime * Year 
##                                             Value Std.Error  DF    t-value
## Thermal_regimeTropical upwelling        -325.2603  416.3072   7 -0.7812988
## Thermal_regimeEquatorial upwelling       233.0673  642.0373   7  0.3630121
## Thermal_regimeSeasonal                  1091.9919  953.7136   7  1.1449893
## Year                                       0.1794    0.2081 116  0.8620877
## Thermal_regimeEquatorial upwelling:Year   -0.2869    0.3827 116 -0.7494770
## Thermal_regimeSeasonal:Year               -0.7108    0.5204 116 -1.3659955
##                                         p-value
## Thermal_regimeTropical upwelling         0.4602
## Thermal_regimeEquatorial upwelling       0.7273
## Thermal_regimeSeasonal                   0.2898
## Year                                     0.3904
## Thermal_regimeEquatorial upwelling:Year  0.4551
## Thermal_regimeSeasonal:Year              0.1746
##  Correlation: 
##                                         Thr_Tu Thr_Eu Thrm_S Year   T_Eu:Y
## Thermal_regimeEquatorial upwelling       0.000                            
## Thermal_regimeSeasonal                   0.000  0.000                     
## Year                                    -1.000  0.000  0.000              
## Thermal_regimeEquatorial upwelling:Year  0.544 -0.839  0.000 -0.544       
## Thermal_regimeSeasonal:Year              0.400  0.000 -0.917 -0.400  0.217
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -2.0543416 -0.6676586 -0.0749493  0.6290682  2.7690585 
## 
## Number of Observations: 128
## Number of Groups: 10
  anova(model1_Model1_NoThermallyStable)
Model1_NoThermallyStable <- ggplot(data=subset(aggr.region, Thermal_regime!="Thermally stable"), aes(x=Year, y=Coral_cover)) +
  geom_boxplot(aes(group=Year), outlier.shape = NA) +
  stat_boxplot(aes(group=Year), geom = 'errorbar')+
  geom_smooth(method = lm, se=FALSE, linetype = "dashed", colour="red")+
  geom_smooth(span = 0.2, se=T, colour="darkgray")+ 
  #geom_point(aes(colour=Thermal_regime), alpha=0.5) +
  scale_y_continuous("Live coral cover (%)", limits = c(-5, 90),
                     breaks = seq(0, 90, by=20), expand = c(0,0))+
  scale_x_continuous("", limits = c(1969, 2015),
                     breaks = seq(1970, 2014, by=2), expand = c(0,0))+
  scale_color_manual(values=my_colours) + labs(title="a")+
  #annotate("text", x = 1983, y = 70, label = "*", size=8)+
  annotate("rect", xmin = 1982, xmax = 1983, ymin = 0, ymax = 80,  alpha = .2, fill="red")+
  annotate("rect", xmin = 1997, xmax = 1998, ymin = 0, ymax = 80,  alpha = .2, fill="red")+
  annotate("rect", xmin = 1972, xmax = 1973, ymin = 0, ymax = 80,  alpha = .2, fill="red")
# Model1_NoThermallyStable
  • Remove tropical upwelling
 model1_Model1_NoTropical_Upwelling <- lme(
    Coral_cover ~ -1 + Thermal_regime * Year, random = ~1|Region,
    data=subset(aggr.region, Thermal_regime!="Tropical upwelling"))
  summary(model1_Model1_NoTropical_Upwelling)
## Linear mixed-effects model fit by REML
##  Data: subset(aggr.region, Thermal_regime != "Tropical upwelling") 
##        AIC      BIC    logLik
##   1071.526 1093.958 -527.7628
## 
## Random effects:
##  Formula: ~1 | Region
##         (Intercept) Residual
## StdDev:    7.304633 15.29287
## 
## Fixed effects: Coral_cover ~ -1 + Thermal_regime * Year 
##                                             Value Std.Error  DF    t-value
## Thermal_regimeThermally stable          -383.7268  342.7238 112 -1.1196384
## Thermal_regimeEquatorial upwelling       186.0471  526.2945  11  0.3535038
## Thermal_regimeSeasonal                   996.7332  763.3534 112  1.3057296
## Year                                       0.2043    0.1715 112  1.1908313
## Thermal_regimeEquatorial upwelling:Year   -0.2878    0.3144 112 -0.9156186
## Thermal_regimeSeasonal:Year               -0.6872    0.4157 112 -1.6532810
##                                         p-value
## Thermal_regimeThermally stable           0.2653
## Thermal_regimeEquatorial upwelling       0.7304
## Thermal_regimeSeasonal                   0.1943
## Year                                     0.2362
## Thermal_regimeEquatorial upwelling:Year  0.3618
## Thermal_regimeSeasonal:Year              0.1011
##  Correlation: 
##                                         Thr_Ts Thr_Eu Thrm_S Year   T_Eu:Y
## Thermal_regimeEquatorial upwelling       0.000                            
## Thermal_regimeSeasonal                   0.016  0.000                     
## Year                                    -1.000  0.000 -0.017              
## Thermal_regimeEquatorial upwelling:Year  0.546 -0.838  0.009 -0.546       
## Thermal_regimeSeasonal:Year              0.398  0.000 -0.911 -0.397  0.217
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.17613688 -0.78949307 -0.01633052  0.53848693  2.59536164 
## 
## Number of Observations: 128
## Number of Groups: 12
  anova(model1_Model1_NoTropical_Upwelling)
Model1_NoTropical_Upwelling <- ggplot(data=subset(aggr.region, Thermal_regime!="Tropical upwelling"), aes(x=Year, y=Coral_cover)) +
  geom_boxplot(aes(group=Year), outlier.shape = NA) +
  stat_boxplot(aes(group=Year), geom = 'errorbar')+
  geom_smooth(method = lm, se=FALSE, linetype = "dashed", colour="red")+
  geom_smooth(span = 0.2, se=T, colour="darkgray")+ 
  #geom_point(aes(colour=Thermal_regime), alpha=0.5) +
  scale_y_continuous("Live coral cover (%)", limits = c(-5, 90), 
                     breaks = seq(0, 90, by=20), expand = c(0,0))+
  scale_x_continuous("", limits = c(1969, 2015),
                     breaks = seq(1970, 2014, by=2), expand = c(0,0))+
  scale_color_manual(values=my_colours) + labs(title="b")+
  #annotate("text", x = 1983, y = 70, label = "*", size=8)+
  annotate("rect", xmin = 1982, xmax = 1983, ymin = 0, ymax = 80,  alpha = .2, fill="red")+
  annotate("rect", xmin = 1997, xmax = 1998, ymin = 0, ymax = 80,  alpha = .2, fill="red")+
  annotate("rect", xmin = 1972, xmax = 1973, ymin = 0, ymax = 80,  alpha = .2, fill="red")
#Model1_NoTropical_Upwelling
  • Remove equatorial upwelling
 model1_Model1_NoEquatorialUpwelling <- lme(
    Coral_cover ~ -1 + Thermal_regime * Year, random = ~1|Region,
    data=subset(aggr.region, Thermal_regime!="Equatorial upwelling"))
  summary(model1_Model1_NoEquatorialUpwelling)
## Linear mixed-effects model fit by REML
##  Data: subset(aggr.region, Thermal_regime != "Equatorial upwelling") 
##      AIC     BIC    logLik
##   1519.4 1544.58 -751.7001
## 
## Random effects:
##  Formula: ~1 | Region
##         (Intercept) Residual
## StdDev:    13.05005 16.25116
## 
## Fixed effects: Coral_cover ~ -1 + Thermal_regime * Year 
##                                           Value Std.Error  DF    t-value
## Thermal_regimeThermally stable        -377.3963  372.0442 161 -1.0143858
## Thermal_regimeTropical upwelling      -331.3816  366.4117  12 -0.9043969
## Thermal_regimeSeasonal                1044.5202  831.3348 161  1.2564375
## Year                                     0.2009    0.1861 161  1.0797967
## Thermal_regimeTropical upwelling:Year   -0.0185    0.2611 161 -0.0710180
## Thermal_regimeSeasonal:Year             -0.7077    0.4527 161 -1.5632421
##                                       p-value
## Thermal_regimeThermally stable         0.3119
## Thermal_regimeTropical upwelling       0.3836
## Thermal_regimeSeasonal                 0.2108
## Year                                   0.2818
## Thermal_regimeTropical upwelling:Year  0.9435
## Thermal_regimeSeasonal:Year            0.1200
##  Correlation: 
##                                       Thr_Ts Thr_Tu Thrm_S Year   T_Tu:Y
## Thermal_regimeTropical upwelling       0.000                            
## Thermal_regimeSeasonal                 0.013  0.000                     
## Year                                  -1.000  0.000 -0.014              
## Thermal_regimeTropical upwelling:Year  0.713 -0.701  0.010 -0.713       
## Thermal_regimeSeasonal:Year            0.399  0.000 -0.912 -0.398  0.283
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.33867507 -0.65989684 -0.01749346  0.63243168  3.15005831 
## 
## Number of Observations: 178
## Number of Groups: 13
  anova(model1_Model1_NoEquatorialUpwelling)
Model1_No_EqUpwelling <- ggplot(data=subset(aggr.region, Thermal_regime!="Equatorial upwelling"), aes(x=Year, y=Coral_cover)) +
  geom_boxplot(aes(group=Year), outlier.shape = NA) +
  stat_boxplot(aes(group=Year), geom = 'errorbar')+
  geom_smooth(method = lm, se=FALSE, linetype = "dashed", colour="red")+
  geom_smooth(span = 0.2, se=T, colour="darkgray")+ 
  #geom_point(aes(colour=Thermal_regime), alpha=0.5) +
  scale_y_continuous("Live coral cover (%)", limits = c(-5, 90), 
                     breaks = seq(0, 90, by=20), expand = c(0,0))+
  scale_x_continuous("", limits = c(1969, 2015),
                     breaks = seq(1970, 2014, by=2), expand = c(0,0))+
  scale_color_manual(values=my_colours) + labs(title="c")+
  #annotate("text", x = 1983, y = 70, label = "*", size=8)+
  annotate("rect", xmin = 1982, xmax = 1983, ymin = 0, ymax = 80,  alpha = .2, fill="red")+
  annotate("rect", xmin = 1997, xmax = 1998, ymin = 0, ymax = 80,  alpha = .2, fill="red")+
  annotate("rect", xmin = 1972, xmax = 1973, ymin = 0, ymax = 80,  alpha = .2, fill="red")
# Model1_No_EqUpwelling
  • Remove seasonal
 model1_Model1_NoSeasonal <- lme(
    Coral_cover ~ -1 + Thermal_regime * Year, random = ~1|Region,
    data=subset(aggr.region, Thermal_regime!="Seasonal"))
  summary(model1_Model1_NoSeasonal)
## Linear mixed-effects model fit by REML
##  Data: subset(aggr.region, Thermal_regime != "Seasonal") 
##        AIC      BIC    logLik
##   1501.822 1526.861 -742.9109
## 
## Random effects:
##  Formula: ~1 | Region
##         (Intercept) Residual
## StdDev:    13.68379  16.5809
## 
## Fixed effects: Coral_cover ~ -1 + Thermal_regime * Year 
##                                             Value Std.Error  DF    t-value
## Thermal_regimeThermally stable          -386.8361  380.4004  11 -1.0169181
## Thermal_regimeTropical upwelling        -336.7578  374.1751  11 -0.9000007
## Thermal_regimeEquatorial upwelling       236.1062  576.4625  11  0.4095777
## Year                                       0.2060    0.1903 159  1.0825575
## Thermal_regimeTropical upwelling:Year     -0.0210    0.2668 159 -0.0786044
## Thermal_regimeEquatorial upwelling:Year   -0.3151    0.3455 159 -0.9118227
##                                         p-value
## Thermal_regimeThermally stable           0.3310
## Thermal_regimeTropical upwelling         0.3874
## Thermal_regimeEquatorial upwelling       0.6900
## Year                                     0.2806
## Thermal_regimeTropical upwelling:Year    0.9374
## Thermal_regimeEquatorial upwelling:Year  0.3632
##  Correlation: 
##                                         Thr_Ts Thr_Tu Thr_Eu Year   T_Tu:Y
## Thermal_regimeTropical upwelling         0.000                            
## Thermal_regimeEquatorial upwelling       0.000  0.000                     
## Year                                    -1.000  0.000  0.000              
## Thermal_regimeTropical upwelling:Year    0.713 -0.701  0.000 -0.713       
## Thermal_regimeEquatorial upwelling:Year  0.551  0.000 -0.835 -0.551  0.393
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.29410007 -0.67127700 -0.01898876  0.68352406  3.08803534 
## 
## Number of Observations: 175
## Number of Groups: 14
  anova(model1_Model1_NoSeasonal)
Model1_NoSeasonal <- ggplot(data=subset(aggr.region, Thermal_regime!="Seasonal"), aes(x=Year, y=Coral_cover)) +
  geom_boxplot(aes(group=Year), outlier.shape = NA) +
  stat_boxplot(aes(group=Year), geom = 'errorbar')+
  geom_smooth(method = lm, se=FALSE, linetype = "dashed", colour="red")+
  geom_smooth(span = 0.2, se=T, colour="darkgray")+ 
  #geom_point(aes(colour=Thermal_regime), alpha=0.5) +
  scale_y_continuous("Live coral cover (%)", limits = c(-5, 90),
                     breaks = seq(0, 90, by=20), expand = c(0,0))+
  scale_x_continuous("", limits = c(1969, 2015),
                     breaks = seq(1970, 2014, by=2), expand = c(0,0))+
  scale_color_manual(values=my_colours) + labs(title="d")+
  #annotate("text", x = 1983, y = 70, label = "*", size=8)+
  annotate("rect", xmin = 1982, xmax = 1983, ymin = 0, ymax = 80,  alpha = .2, fill="red")+
  annotate("rect", xmin = 1997, xmax = 1998, ymin = 0, ymax = 80,  alpha = .2, fill="red")+
  annotate("rect", xmin = 1972, xmax = 1973, ymin = 0, ymax = 80,  alpha = .2, fill="red")
#Model1_NoSeasonal

FIGURE S4: Effect of thermal regime removal/addition

FigureS4<-grid.arrange(Model1_NoThermallyStable, Model1_NoTropical_Upwelling,
                       Model1_No_EqUpwelling, Model1_NoSeasonal, 
                       ncol = 2)

#ggsave(file="FigS4.pdf", plot=FigureS5, width=6.69, height=6.69)

Figure S4: Coral cover trends excluding (a) Thermally stable, (b) Tropical upwelling, (c) Equatorial upwelling, and (d) Seasonal thermal regimes

Add CI based on bootstrapping to fig 2?

To use “bootMer” models need to be run with lme4 instead of nlme. Please check accuracy of the new (b) models

Model 1b: All years 1970_2014 with lme4

  # 1. Build lme4 model:
   
    model1b <- lmer(
      Coral_cover ~ Thermal_regime * Year + 
        (1|Region), data=aggr.region, REML = FALSE)
  
    summary(model1b)
## Linear mixed model fit by maximum likelihood . t-tests use
##   Satterthwaite's method [lmerModLmerTest]
## Formula: Coral_cover ~ Thermal_regime * Year + (1 | Region)
##    Data: aggr.region
## 
##      AIC      BIC   logLik deviance df.resid 
##   1757.4   1790.5   -868.7   1737.4      193 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.27927 -0.72367 -0.01805  0.63561  3.09102 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Region   (Intercept)  95.09    9.751  
##  Residual             272.58   16.510  
## Number of obs: 203, groups:  Region, 16
## 
## Fixed effects:
##                                           Estimate Std. Error         df
## (Intercept)                             -376.66123  373.68559  201.99026
## Thermal_regimeTropical upwelling         113.26629  524.31829  202.24674
## Thermal_regimeEquatorial upwelling       588.01224  682.19129  197.01240
## Thermal_regimeSeasonal                  1389.07741  907.62131  202.83558
## Year                                       0.20069    0.18698  201.75471
## Thermal_regimeTropical upwelling:Year     -0.05175    0.26224  202.06878
## Thermal_regimeEquatorial upwelling:Year   -0.29705    0.34139  196.60807
## Thermal_regimeSeasonal:Year               -0.69147    0.45351  202.75097
##                                         t value Pr(>|t|)
## (Intercept)                              -1.008    0.315
## Thermal_regimeTropical upwelling          0.216    0.829
## Thermal_regimeEquatorial upwelling        0.862    0.390
## Thermal_regimeSeasonal                    1.530    0.127
## Year                                      1.073    0.284
## Thermal_regimeTropical upwelling:Year    -0.197    0.844
## Thermal_regimeEquatorial upwelling:Year  -0.870    0.385
## Thermal_regimeSeasonal:Year              -1.525    0.129
## 
## Correlation of Fixed Effects:
##             (Intr) Thr_Tu Thr_Eu Thrm_S Year   T_Tu:Y T_Eu:Y
## Thrml_rgmTu -0.713                                          
## Thrml_rgmEu -0.548  0.390                                   
## Thrml_rgmSs -0.397  0.283  0.218                            
## Year        -1.000  0.713  0.548  0.397                     
## Thrml_rTu:Y  0.713 -1.000 -0.391 -0.283 -0.713              
## Thrml_rEu:Y  0.548 -0.390 -1.000 -0.217 -0.548  0.391       
## Thrml_rgS:Y  0.398 -0.283 -0.218 -1.000 -0.397  0.283  0.217
    anova(model1b) # Thermal regime is not significant anymore?? 
    ranova(model1b) # Region seems the only significant effect
    drop1(model1b)

Not sure if this if right since the model itself is not significant 

  # 2. Predict values:
    pred1 <- predict(model1b, re.form = NA)

  #3. Bootstrap CI:
    boot1 <- bootMer(model1b, predict, nsim = 1000, re.form = NULL) # include random effects, reduce CI a lot!
    std.err <- apply(boot1$t, 2, sd)
    CI.lo_1 <- pred1 - std.err*1.96
    CI.hi_1 <- pred1 + std.err*1.96

  #Plot
  Model1b_plot<- ggplot(
    aggr.region, aes(x = Year, y = Coral_cover, colour = Thermal_regime)) +
    geom_line(aes(y = pred1),size=2) +
    geom_point(aes(fill=factor(Thermal_regime)),
             shape = 21, colour = "black", size = 2, stroke = 0.3, alpha=0.5) +
    geom_ribbon(aes(ymin = CI.lo_1, ymax = CI.hi_1),
                size=2, alpha = 0.03, linetype = 0) +
    scale_color_manual(values=my_colours) +
    scale_fill_manual(values=my_colours) +
    scale_y_continuous("Live coral cover (%)", 
                       limits = c(-20, 90), 
                       breaks = seq(0, 90, by=20), expand = c(0,0))+
    scale_x_continuous("", limits = c(1969, 2015),
                     breaks = seq(1970, 2014, by=2), expand = c(0,0))
  
  Model1b_plot

Model 2b: Time-interval 1 - ENSO.73_82 with lme4

Not enough data

Model 3b: Time-interval 2 - ENSO.83_97 with lme4

  # 1. Build lme4 model:
   
    model3b <- lmer(
      Coral_cover ~ Thermal_regime * Year + 
        (1|Region), data=ENSO.83_97, REML = FALSE)
  
    summary(model3b)
## Linear mixed model fit by maximum likelihood . t-tests use
##   Satterthwaite's method [lmerModLmerTest]
## Formula: Coral_cover ~ Thermal_regime * Year + (1 | Region)
##    Data: ENSO.83_97
## 
##      AIC      BIC   logLik deviance df.resid 
##    506.7    528.6   -243.4    486.7       56 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.0485 -0.5244 -0.0721  0.4890  3.2655 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Region   (Intercept) 13.12    3.622   
##  Residual             84.48    9.191   
## Number of obs: 66, groups:  Region, 10
## 
## Fixed effects:
##                                           Estimate Std. Error         df
## (Intercept)                             -2085.2528   849.2084    60.0798
## Thermal_regimeTropical upwelling        -5020.5691  1187.3127    62.4593
## Thermal_regimeEquatorial upwelling       2197.0397  2802.4619    57.5480
## Thermal_regimeSeasonal                   8274.4993  2002.6425    63.4562
## Year                                        1.0546     0.4267    60.0267
## Thermal_regimeTropical upwelling:Year       2.5355     0.5963    62.4179
## Thermal_regimeEquatorial upwelling:Year    -1.1105     1.4097    57.5416
## Thermal_regimeSeasonal:Year                -4.1450     1.0053    63.4138
##                                         t value Pr(>|t|)    
## (Intercept)                              -2.456 0.016978 *  
## Thermal_regimeTropical upwelling         -4.229 7.82e-05 ***
## Thermal_regimeEquatorial upwelling        0.784 0.436274    
## Thermal_regimeSeasonal                    4.132 0.000107 ***
## Year                                      2.472 0.016307 *  
## Thermal_regimeTropical upwelling:Year     4.252 7.22e-05 ***
## Thermal_regimeEquatorial upwelling:Year  -0.788 0.434073    
## Thermal_regimeSeasonal:Year              -4.123 0.000111 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Thr_Tu Thr_Eu Thrm_S Year   T_Tu:Y T_Eu:Y
## Thrml_rgmTu -0.715                                          
## Thrml_rgmEu -0.303  0.217                                   
## Thrml_rgmSs -0.414  0.296  0.125                            
## Year        -1.000  0.715  0.303  0.414                     
## Thrml_rTu:Y  0.716 -1.000 -0.217 -0.296 -0.716              
## Thrml_rEu:Y  0.303 -0.216 -1.000 -0.125 -0.303  0.217       
## Thrml_rgS:Y  0.414 -0.296 -0.126 -1.000 -0.414  0.296  0.125
    anova(model3b) # Thermal regime and its interaction with year are significant 
    ranova(model3b) # random effects are not 
    drop1(model3b)
  # 2. Predict values:
    pred3 <- predict(model3b, re.form = NA)

  #3. Bootstrap CI:
    boot3 <- bootMer(model3b, predict, nsim = 1000, re.form = NA) # No random effects
    std.err <- apply(boot3$t, 2, sd)
    CI.lo_3 <- pred3 - std.err*1.96
    CI.hi_3 <- pred3 + std.err*1.96

  #Plot
    
   Model3b_plot<- ggplot(
     ENSO.83_97, aes(x = Year, y = Coral_cover, colour = Thermal_regime))+
    geom_line(aes(y = pred3),size=2) +
    geom_point(aes(fill=factor(Thermal_regime)),
          shape = 21, colour = "black", size = 2, stroke = 0.3, alpha=0.5) + 
    geom_ribbon(aes(ymin = CI.lo_3, ymax = CI.hi_3),
                size=2, alpha = 0.03, linetype = 0) +
    scale_fill_manual(values=my_colours)+
    scale_colour_manual(values=my_colours) +
    scale_x_continuous("", limits = c(1982.7, 1998.2), 
                        breaks = seq(1983, 1997, by=2), 
                        expand = c(0.02,0.02))+
    scale_y_continuous("Live coral cover (%)", 
                       limits = c(-20, 90), 
                       breaks = seq(0, 90, by=20), expand = c(0,0))
    
   Model3b_plot

  #ggsave(file="Model3b_plot.png", plot=Model3b_plot, width=4, height=4)
  #ggsave(file="Model3b_plot.pdf", plot=Model3b_plot, width=4, height=4)

Model 4b: Time-interval 3 - ENSO.98_14 with lme4

  # 1. Build lme4 model:
   
    model4b <- lmer(
      Coral_cover ~ Thermal_regime * Year + 
        (1|Region), data=ENSO.98_14, REML = FALSE)
  
    summary(model4b)
## Linear mixed model fit by maximum likelihood . t-tests use
##   Satterthwaite's method [lmerModLmerTest]
## Formula: Coral_cover ~ Thermal_regime * Year + (1 | Region)
##    Data: ENSO.98_14
## 
##      AIC      BIC   logLik deviance df.resid 
##    896.3    922.8   -438.2    876.3       94 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.9576 -0.5586 -0.0027  0.2837  2.7825 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Region   (Intercept)  81.54    9.03   
##  Residual             228.45   15.11   
## Number of obs: 104, groups:  Region, 12
## 
## Fixed effects:
##                                           Estimate Std. Error         df
## (Intercept)                             -1723.2489  1349.0589   103.8734
## Thermal_regimeTropical upwelling         4931.1394  1694.6395   103.9964
## Thermal_regimeEquatorial upwelling       -397.0246  3139.2328    96.2899
## Thermal_regimeSeasonal                   1958.6900  2050.7623   100.1819
## Year                                        0.8727     0.6728   103.8861
## Thermal_regimeTropical upwelling:Year      -2.4541     0.8449   103.9980
## Thermal_regimeEquatorial upwelling:Year     0.1964     1.5650    96.2689
## Thermal_regimeSeasonal:Year                -0.9740     1.0221   100.0849
##                                         t value Pr(>|t|)   
## (Intercept)                              -1.277  0.20432   
## Thermal_regimeTropical upwelling          2.910  0.00442 **
## Thermal_regimeEquatorial upwelling       -0.126  0.89962   
## Thermal_regimeSeasonal                    0.955  0.34182   
## Year                                      1.297  0.19747   
## Thermal_regimeTropical upwelling:Year    -2.904  0.00449 **
## Thermal_regimeEquatorial upwelling:Year   0.126  0.90037   
## Thermal_regimeSeasonal:Year              -0.953  0.34293   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Thr_Tu Thr_Eu Thrm_S Year   T_Tu:Y T_Eu:Y
## Thrml_rgmTu -0.796                                          
## Thrml_rgmEu -0.430  0.342                                   
## Thrml_rgmSs -0.651  0.519  0.280                            
## Year        -1.000  0.796  0.430  0.651                     
## Thrml_rTu:Y  0.796 -1.000 -0.342 -0.519 -0.796              
## Thrml_rEu:Y  0.430 -0.342 -1.000 -0.280 -0.430  0.342       
## Thrml_rgS:Y  0.651 -0.519 -0.280 -1.000 -0.651  0.519  0.280
    anova(model4b) # Thermal regime and its interaction with year are significant 
    ranova(model4b) # random effects are significant
    drop1(model4b)
  # 2. Predict values:
    pred4 <- predict(model4b, re.form = NA)

  #3. Bootstrap CI:
    boot4 <- bootMer(model4b, predict, nsim = 1000, re.form = NULL) # Include random effects
    std.err <- apply(boot4$t, 2, sd)
    CI.lo_4 <- pred4 - std.err*1.96
    CI.hi_4 <- pred4 + std.err*1.96

  #Plot
   Model4b_plot<- ggplot(ENSO.98_14, aes(
     x = Year, y = Coral_cover, colour = Thermal_regime)) +
    geom_line(aes(y = pred4),size=2) +
    geom_point(aes(fill=factor(Thermal_regime)),
             shape = 21, colour = "black", size = 2, 
             stroke = 0.3, alpha=0.5) +
    geom_ribbon(aes(ymin = CI.lo_4, ymax = CI.hi_4),size=2, 
                alpha = 0.03, linetype = 0) +
    scale_y_continuous("Live coral cover (%)",
                       limits = c(-20, 90), 
                       breaks = seq(0, 90, by=20), expand = c(0,0))+
    scale_x_continuous("", limits = c(1997.7, 2014.3),
                     breaks = seq(1998, 2014, by=2), expand = c(0.02,0.02))+
   scale_color_manual(values=my_colours) +
   scale_fill_manual(values=my_colours)
   Model4b_plot

  #ggsave(file="Model4b_plot.png", plot=Model4b_plot, width=4, height=4)
  #ggsave(file="Model4b_plot.pdf", plot=Model4b_plot, width=4, height=4)

Model 5b: Time-interval 4 - ENSO.98_14 with lme4

  # 1. Build lme4 model:
   
    model5b <- lmer(
      Coral_cover ~ Thermal_regime * Year + 
        (1|Region), data=ENSO.83_14, REML = FALSE)
  
    summary(model5b)
## Linear mixed model fit by maximum likelihood . t-tests use
##   Satterthwaite's method [lmerModLmerTest]
## Formula: Coral_cover ~ Thermal_regime * Year + (1 | Region)
##    Data: ENSO.83_14
## 
##      AIC      BIC   logLik deviance df.resid 
##   1453.3   1484.7   -716.7   1433.3      160 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.4879 -0.5122 -0.1008  0.4671  3.0884 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Region   (Intercept)  95.58    9.777  
##  Residual             235.87   15.358  
## Number of obs: 170, groups:  Region, 13
## 
## Fixed effects:
##                                           Estimate Std. Error         df
## (Intercept)                             -1793.4643   507.5993   163.9918
## Thermal_regimeTropical upwelling         1908.5107   710.6627   168.1130
## Thermal_regimeEquatorial upwelling        334.0786  1027.8011   166.7244
## Thermal_regimeSeasonal                   2795.1759   921.5719   168.5964
## Year                                        0.9080     0.2538   163.7270
## Thermal_regimeTropical upwelling:Year      -0.9480     0.3551   167.9676
## Thermal_regimeEquatorial upwelling:Year    -0.1683     0.5134   166.5944
## Thermal_regimeSeasonal:Year                -1.3932     0.4604   168.3969
##                                         t value Pr(>|t|)    
## (Intercept)                              -3.533 0.000533 ***
## Thermal_regimeTropical upwelling          2.686 0.007967 ** 
## Thermal_regimeEquatorial upwelling        0.325 0.745557    
## Thermal_regimeSeasonal                    3.033 0.002804 ** 
## Year                                      3.578 0.000455 ***
## Thermal_regimeTropical upwelling:Year    -2.670 0.008338 ** 
## Thermal_regimeEquatorial upwelling:Year  -0.328 0.743526    
## Thermal_regimeSeasonal:Year              -3.026 0.002866 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Thr_Tu Thr_Eu Thrm_S Year   T_Tu:Y T_Eu:Y
## Thrml_rgmTu -0.714                                          
## Thrml_rgmEu -0.494  0.353                                   
## Thrml_rgmSs -0.536  0.383  0.265                            
## Year        -1.000  0.714  0.494  0.536                     
## Thrml_rTu:Y  0.715 -1.000 -0.353 -0.383 -0.715              
## Thrml_rEu:Y  0.494 -0.353 -1.000 -0.265 -0.494  0.353       
## Thrml_rgS:Y  0.536 -0.383 -0.265 -1.000 -0.536  0.383  0.265
    anova(model5b) # Thermal regime and its interaction with year are significant 
    ranova(model5b) # random effects are significant
    drop1(model5b)
  # 2. Predict values:
    pred5 <- predict(model5b,re.form = NA)

  #3. Bootstrap CI:
    boot5 <- bootMer(model5b, predict, nsim = 1000, re.form = NULL) # Include random effects
    std.err <- apply(boot5$t, 2, sd)
    CI.lo_5 <- pred5 - std.err*1.96
    CI.hi_5 <- pred5 + std.err*1.96

  #Plot
   Model5b_plot<- ggplot(ENSO.83_14, aes(
     x = Year, y = Coral_cover, colour = Thermal_regime)) +
    geom_line(aes(y = pred5),size=2) +
    geom_point(aes(fill=factor(Thermal_regime)),
             shape = 21, colour = "black", size = 2, 
             stroke = 0.3, alpha=0.5) +
    geom_ribbon(aes(ymin = CI.lo_5, ymax = CI.hi_5),size=2, 
                alpha = 0.03, linetype = 0) +
    scale_y_continuous("Live coral cover (%)",
                       limits = c(-20, 90), 
                       breaks = seq(0, 90, by=20), expand = c(0,0))+
    scale_x_continuous("", limits = c(1983, 2014.3),
                     breaks = seq(1983, 2014, by=2), expand = c(0.02,0.02))+
   scale_color_manual(values=my_colours) +
   scale_fill_manual(values=my_colours)
   Model5b_plot

  #ggsave(file="Model5b_plot.png", plot=Model4b_plot, width=4, height=4)

5. Packages used

# Creates bibliography 
#knitr::write_bib(c(.packages()), "packages.bib")

Auguie, Baptiste. 2017. GridExtra: Miscellaneous Functions for “Grid” Graphics. https://CRAN.R-project.org/package=gridExtra.

Bartoń, Kamil. 2019. MuMIn: Multi-Model Inference. https://CRAN.R-project.org/package=MuMIn.

Bates, Douglas, and Martin Maechler. 2019. Matrix: Sparse and Dense Matrix Classes and Methods. https://CRAN.R-project.org/package=Matrix.

Bates, Douglas, Martin Maechler, Ben Bolker, and Steven Walker. 2019. Lme4: Linear Mixed-Effects Models Using ’Eigen’ and S4. https://CRAN.R-project.org/package=lme4.

Dowle, Matt, and Arun Srinivasan. 2019. Data.table: Extension of ‘Data.frame‘. https://CRAN.R-project.org/package=data.table.

Edwards, Stefan McKinnon. 2019. Lemon: Freshing up Your ’Ggplot2’ Plots. https://CRAN.R-project.org/package=lemon.

Genz, Alan, Frank Bretz, Tetsuhisa Miwa, Xuefei Mi, and Torsten Hothorn. 2019. Mvtnorm: Multivariate Normal and T Distributions. https://CRAN.R-project.org/package=mvtnorm.

Henry, Lionel, and Hadley Wickham. 2019. Purrr: Functional Programming Tools. https://CRAN.R-project.org/package=purrr.

Hothorn, Torsten. 2019. TH.data: TH’s Data Archive. https://CRAN.R-project.org/package=TH.data.

Hothorn, Torsten, Frank Bretz, and Peter Westfall. 2019. Multcomp: Simultaneous Inference in General Parametric Models. https://CRAN.R-project.org/package=multcomp.

Kuznetsova, Alexandra, Per Bruun Brockhoff, and Rune Haubo Bojesen Christensen. 2019. LmerTest: Tests in Linear Mixed Effects Models. https://CRAN.R-project.org/package=lmerTest.

Müller, Kirill, and Hadley Wickham. 2019. Tibble: Simple Data Frames. https://CRAN.R-project.org/package=tibble.

Neuwirth, Erich. 2014. RColorBrewer: ColorBrewer Palettes. https://CRAN.R-project.org/package=RColorBrewer.

Pinheiro, José, Douglas Bates, and R-core. 2019. Nlme: Linear and Nonlinear Mixed Effects Models. https://CRAN.R-project.org/package=nlme.

R Core Team. 2020. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.

Richard A. Becker, Original S code by, Allan R. Wilks. R version by Ray Brownrigg. Enhancements by Thomas P Minka, and Alex Deckmyn. 2018. Maps: Draw Geographical Maps. https://CRAN.R-project.org/package=maps.

Ripley, Brian. 2019. MASS: Support Functions and Datasets for Venables and Ripley’s Mass. https://CRAN.R-project.org/package=MASS.

Therneau, Terry M. 2019. Survival: Survival Analysis. https://CRAN.R-project.org/package=survival.

Wickham, Hadley. 2017a. Reshape2: Flexibly Reshape Data: A Reboot of the Reshape Package. https://CRAN.R-project.org/package=reshape2.

———. 2017b. Tidyverse: Easily Install and Load the ’Tidyverse’. https://CRAN.R-project.org/package=tidyverse.

———. 2019a. Forcats: Tools for Working with Categorical Variables (Factors). https://CRAN.R-project.org/package=forcats.

———. 2019b. Stringr: Simple, Consistent Wrappers for Common String Operations. https://CRAN.R-project.org/package=stringr.

Wickham, Hadley, and Lionel Henry. 2020. Tidyr: Tidy Messy Data. https://CRAN.R-project.org/package=tidyr.

Wickham, Hadley, Winston Chang, Lionel Henry, Thomas Lin Pedersen, Kohske Takahashi, Claus Wilke, Kara Woo, and Hiroaki Yutani. 2019. Ggplot2: Create Elegant Data Visualisations Using the Grammar of Graphics. https://CRAN.R-project.org/package=ggplot2.

Wickham, Hadley, Romain François, Lionel Henry, and Kirill Müller. 2019. Dplyr: A Grammar of Data Manipulation. https://CRAN.R-project.org/package=dplyr.

Wickham, Hadley, Jim Hester, and Romain Francois. 2018. Readr: Read Rectangular Text Data. https://CRAN.R-project.org/package=readr.

Wright, Kevin. 2018. Pals: Color Palettes, Colormaps, and Tools to Evaluate Them. https://CRAN.R-project.org/package=pals.