Uva Island bleaching v2
Ana M. Palacio-Castro and Tyler Smith
2022-08-30
Setup
1. Binary approach - UVA ISLAND
Data
Counts1983_1997<-read.csv("BleachingData/Counts_1983_1997.csv",
header = T)
# summary(Counts1983_1997)
# sum(Counts1983_1997$Number)
# sum(table(rep(1:nrow(Counts1983_1997), Counts1983_1997$Number)))
Bi_1983_97 <- Counts1983_1997[rep(1:nrow(Counts1983_1997),
Counts1983_1997$Number), 1:4]
#nrow(Bi_1983_97)
Bi_1983_97$Site<-NA
Bi_2014_16<-read.csv("BleachingData/Counts_2014_2016.csv",
header = T)
# nrow(Bi_2014_16)
# summary(Bi_2014_16)
Binomial_all.data<-rbind(Bi_1983_97, Bi_2014_16)
#summary(Binomial_all.data)
#write.csv(Binomial.data, "BleachingData/Binomial_data.csv", row.names = F)#Binomial.data<-read.csv("BleachingData/Binomial_data.csv", header = T)
Species<-read.csv("BleachingData/Genus.csv", header = T)
Binomial_all.data$Date<-as.Date(Binomial_all.data$Date,
format="%m/%d/%y")
Binomial_all.data$Year<-year(Binomial_all.data$Date)
Binomial_all.data$Year_F<-as.factor(Binomial_all.data$Year)
Binomial_all.data<-merge(Binomial_all.data, Species, by="Species", all.x=TRUE)
Binomial_all.data$Sensitivity<-factor(Binomial_all.data$Sensitivity,
levels = c("Millepora",
"Other", "Pocillopora"))
Binomial_all.data$Genus<-factor(Binomial_all.data$Genus,
levels = c("Millepora",
"Pavona", "Gardineroseris",
"Porites",
"Psammocora", "Pocillopora"))
Binomial_all.data$Species[
Binomial_all.data$Species=="Psammocora superficialis"]<-"Psammocora stellata"
Binomial_all.data$Species<-factor(Binomial_all.data$Species, levels = c(
"Millepora platyphylla", "Millepora intricata",
"Pavona chiriquiensis", "Pavona varians", "Pavona clavus", "Porites panamensis",
"Pocillopora elegans", "Pocillopora damicornis", "Pocillopora",
"Gardineroseris planulata",
"Porites lobata", "Pavona gigantea", "Psammocora stellata"))
Binomial_all.data$Affected[Binomial_all.data$Condition=="Healthy"]<-0
Binomial_all.data$Affected[Binomial_all.data$Condition=="Affected"]<-1
# filter data
# Remove 2018 - Recovery
Binomial_all.data<-Binomial_all.data[Binomial_all.data$Year<2018,]
# Remove - species with low replication <20 colonies by year
Binomial_all.data<-Binomial_all.data[Binomial_all.data$Species!="Pavona gigantea",]
Binomial_all.data<-Binomial_all.data[Binomial_all.data$Species!="Millepora platyphylla",]
Binomial_all.data<-Binomial_all.data[Binomial_all.data$Species!="Pavona chiriquiensis",]
Binomial_all.data<-Binomial_all.data[Binomial_all.data$Species!="Psammocora stellata",]
Binomial_all.data<-Binomial_all.data[Binomial_all.data$Species!="Porites panamensis",]
# Remove - non-uva locations
Binomial.data<-Binomial_all.data[Binomial_all.data$Location=="Uva",]
Binomial.data<-(droplevels(Binomial.data))
summary(Binomial.data)## Species Location Date Condition
## Millepora intricata :223 Uva:2294 Min. :1983-03-20 Affected:1303
## Pavona varians :495 1st Qu.:1983-03-20 Healthy : 991
## Pavona clavus :249 Median :1997-10-14
## Pocillopora elegans :291 Mean :1997-02-06
## Pocillopora damicornis :557 3rd Qu.:1998-03-13
## Gardineroseris planulata:351 Max. :2016-04-21
## Porites lobata :128
## Site Year Year_F Genus
## Length:2294 Min. :1983 1983:779 Millepora :223
## Class :character 1st Qu.:1983 1997:558 Pavona :744
## Mode :character Median :1997 1998:421 Gardineroseris:351
## Mean :1997 2014: 80 Porites :128
## 3rd Qu.:1998 2015:271 Pocillopora :848
## Max. :2016 2016:185
##
## Sensitivity Affected
## Millepora : 223 Min. :0.000
## Other :1223 1st Qu.:0.000
## Pocillopora: 848 Median :1.000
## Mean :0.568
## 3rd Qu.:1.000
## Max. :1.000
## Binomial_other<-Binomial_all.data[Binomial_all.data$Location!="Uva",]
Binomial_other<-(droplevels(Binomial_other))
summary(Binomial_other)## Species Location Date
## Millepora intricata : 4 Jicaron-Jicarita:760 Min. :1997-10-16
## Pavona varians : 60 Silva de Afuera :220 1st Qu.:1998-03-11
## Pavona clavus :234 Canal de Afuera : 44 Median :1998-03-11
## Pocillopora elegans :221 Secas :168 Mean :2001-04-09
## Pocillopora damicornis :175 3rd Qu.:1998-03-11
## Gardineroseris planulata:244 Max. :2016-04-20
## Porites lobata :254
## Condition Site Year Year_F
## Affected:332 Length:1192 Min. :1997 1997:220
## Healthy :860 Class :character 1st Qu.:1998 1998:760
## Mode :character Median :1998 2014: 1
## Mean :2001 2015:117
## 3rd Qu.:1998 2016: 94
## Max. :2016
##
## Genus Sensitivity Affected
## Millepora : 4 Millepora : 4 Min. :0.0000
## Pavona :294 Other :792 1st Qu.:0.0000
## Gardineroseris:244 Pocillopora:396 Median :0.0000
## Porites :254 Mean :0.2785
## Pocillopora :396 3rd Qu.:1.0000
## Max. :1.0000
## Number of colonies sampled -Uva
Heat stress at the time of the surveys
DHW_Daily<-read.csv("BleachingData/DHW_Daily_CA.csv", header = T)
summary(DHW_Daily)## Location lat lon X.lon
## Canal de Afuera :12053 Min. :7.242 Min. :277.8 Min. :-82.23
## Jicarón and Jicarita:12053 1st Qu.:7.476 1st Qu.:278.0 1st Qu.:-82.07
## Montuoso :12053 Median :7.793 Median :278.2 Median :-81.91
## Secas :12053 Mean :7.712 Mean :278.1 Mean :-81.96
## Silva de Afuera :12053 3rd Qu.:7.952 3rd Qu.:278.2 3rd Qu.:-81.81
## Uva :12896 Max. :8.090 Max. :278.2 Max. :-81.79
## NA's :12053 NA's :12053 NA's :24949
## SST Date Year MMM
## Min. :24.41 2016-12-19: 7 Min. :1982 Min. :28.94
## 1st Qu.:27.97 2016-12-20: 7 1st Qu.:1991 1st Qu.:29.16
## Median :28.58 2016-12-21: 7 Median :1999 Median :29.16
## Mean :28.55 2016-12-22: 7 Mean :1999 Mean :29.17
## 3rd Qu.:29.17 2016-12-23: 7 3rd Qu.:2008 3rd Qu.:29.25
## Max. :32.00 2016-12-24: 7 Max. :2019 Max. :29.29
## (Other) :73119
## HotSpot D_Stress W_Stress NOAA_DHW_OI
## Min. :0.0000 Min. :0.00000 Min. :0.000000 Min. : 0.0000
## 1st Qu.:0.0000 1st Qu.:0.00000 1st Qu.:0.000000 1st Qu.: 0.0000
## Median :0.0000 Median :0.00000 Median :0.000000 Median : 0.0000
## Mean :0.1207 Mean :0.03748 Mean :0.005354 Mean : 0.4507
## 3rd Qu.:0.0000 3rd Qu.:0.00000 3rd Qu.:0.000000 3rd Qu.: 0.0000
## Max. :2.8355 Max. :2.83555 Max. :0.405078 Max. :12.8800
## NA's :498Location_date<-expand.grid(Location=unique(Binomial.data$Location),
Date=unique(Binomial.data$Date))
Dates<-expand.grid(Date=unique(Binomial.data$Date))
DHW_Daily$Date<-as.Date(DHW_Daily$Date)
DHW_month<-DHW_Daily %>%
group_by(Location) %>%
mutate(DHW_max=rollapplyr(NOAA_DHW_OI,list(-(180:0)),
max,fill=NA, partial=FALSE)) %>% ungroup
#
# DHW_month<-DHW_Daily %>%
# group_by(Location) %>%
# mutate(DHW_max=rollapplyr(NOAA_DHW_OI,list(-(180:0)),
# max,fill=NA, partial=FALSE)) %>% ungroup
DHW_month<-as.data.frame(DHW_month)
DHW_month<- DHW_month %>% dplyr::select("Location", "Date","NOAA_DHW_OI", "DHW_max")
Location_date<-merge(Location_date, DHW_month,
by=c("Location", "Date"), all.x=TRUE)
Binomial.data<-merge(Binomial.data, Location_date,
by=c("Location", "Date"), all.x=TRUE)
str(Binomial.data)## 'data.frame': 2294 obs. of 12 variables:
## $ Location : Factor w/ 1 level "Uva": 1 1 1 1 1 1 1 1 1 1 ...
## $ Date : Date, format: "1983-03-20" "1983-03-20" ...
## $ Species : Factor w/ 7 levels "Millepora intricata",..: 6 6 6 6 6 6 6 6 6 6 ...
## $ Condition : Factor w/ 2 levels "Affected","Healthy": 1 2 1 2 2 2 1 2 2 2 ...
## $ Site : chr NA NA NA NA ...
## $ Year : num 1983 1983 1983 1983 1983 ...
## $ Year_F : Factor w/ 6 levels "1983","1997",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Genus : Factor w/ 5 levels "Millepora","Pavona",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ Sensitivity: Factor w/ 3 levels "Millepora","Other",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ Affected : num 1 0 1 0 0 0 1 0 0 0 ...
## $ NOAA_DHW_OI: num 6.09 6.09 6.09 6.09 6.09 ...
## $ DHW_max : num 6.09 6.09 6.09 6.09 6.09 ...DHW vs Max DHW in 6 months
Dates1<-as.data.frame(
seq.Date(as.Date("1982-08-15"), as.Date("1983-12-15"),
by = "day"))
colnames(Dates1)<-"Date"
Dates1$Event<- "E_1982-83"
Dates1<-merge(Dates1, DHW_month, by=c("Date"), all.x=TRUE )
Dates2<-as.data.frame(
seq.Date(as.Date("1997-07-01"), as.Date("1998-07-31"),
by = "day"))
colnames(Dates2)<-"Date"
Dates2$Event<- "E_1997-1998"
Dates2<-merge(Dates2, DHW_month, by=c("Date"), all.x=TRUE )
Dates3<-as.data.frame(
seq.Date(as.Date("2015-05-01"), as.Date("2016-07-31"),
by = "day"))
colnames(Dates3)<-"Date"
Dates3$Event<- "E_2015-2016"
Dates3<-merge(Dates3, DHW_month, by=c("Date"), all.x=TRUE )
# Date<-rbind(Dates1, Dates2, Dates3)
# Date$Date<-as.Date(Date$Date)DHW_83 <- ggplot(Dates1[Dates1$Location=="Uva", ],
aes(x=Date, y=NOAA_DHW_OI)) +
xlab(" ") +
scale_y_continuous(name="DHW",
limits = c(0,10.5))+
#facet_wrap(Event~., scales = "free_x")+
geom_vline(xintercept = c(Dates$Date), linetype=2)+
geom_line()+
geom_line(aes(x=Date, y=DHW_max), colour="red")
#DHW_83
DHW_98 <- ggplot(Dates2[Dates2$Location=="Uva", ],
aes(x=Date, y=NOAA_DHW_OI)) +
xlab(" ") +
scale_y_continuous(name="DHW",
limits = c(0,10.5))+
#facet_wrap(Event~., scales = "free_x")+
geom_vline(xintercept = c(Dates$Date), linetype=2)+
geom_line()+
geom_line(aes(x=Date, y=DHW_max), colour="red")
#DHW_98
DHW_2016 <- ggplot(Dates3[Dates3$Location=="Uva", ],
aes(x=Date, y=NOAA_DHW_OI)) +
xlab(" ") +
scale_y_continuous(name="DHW",
limits = c(0,10.5))+
#facet_wrap(Event~., scales = "free_x")+
geom_vline(xintercept = c(Dates$Date), linetype=2)+
geom_line()+
geom_line(aes(x=Date, y=DHW_max), colour="red")
#DHW_2016
Max_DHW_plot<-grid.arrange(DHW_83, DHW_98, DHW_2016, ncol=3)
#ggsave(file="Outputs/Max_DHW_plot.svg", plot=Max_DHW_plot, width=5, height=2)Plots by Year
Facet species - Figure 2
binomial_smooth <- function(...) {
geom_smooth(method = "glm", method.args = list(family = "binomial"), ...)
}
Figure_Count<-ggplot(data=Binomial.data[Binomial.data$Year!=2014, ],
aes(x=Year_F, y=Affected)) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5))+
scale_x_discrete(name="Year")+
scale_y_continuous(name="Proportion of colonies affected",
limits = c(0,1.1),
breaks = c(0, 0.25, 0.5, .75, 1))+
geom_abline(slope = 0, intercept = 0.5, colour="grey")+
# geom_jitter(aes(x=Year_F, y=Affected,
# colour=as.factor(Year_F)),
# show.legend = F,
# alpha=0.2, size=0.5, width = 0.2, height = 0.1)+
stat_summary(aes(x=Year_F, y=Affected, fill=Year_F),
shape=21, alpha=0.9,
fun.data = "mean_cl_boot", colour = "black")+
scale_colour_manual(values = c("#FDDBC7", "#67A9CF", "#2166ac",
"#EF8A62", "#B2182B"))+
scale_fill_manual(values = c("#FDDBC7", "#67A9CF", "#2166ac",
"#EF8A62", "#B2182B"))
Figure_Count<-Figure_Count + facet_wrap(Species~., ncol = 4)
Figure_Count
Figure_Count + facet_wrap(Genus~.)
Figure_Count + facet_wrap(Sensitivity~.)
#ggsave(file="Outputs/Figure2d.svg", plot=Figure_Count, width=7, height=5.5)Facet years - Figure S2
Figure_Count_spp2<-ggplot(data=Binomial.data[
Binomial.data$Year!=2014,]) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
axis.title.x = element_blank(),
legend.title = element_blank())+
scale_x_discrete(name="Year")+
scale_y_continuous(name="Proportion of colonies affected",
limits = c(0,1.1),
breaks = c(0, 0.25, 0.5, .75, 1))+
# geom_jitter(aes(x=Species, y=Affected,
# colour=as.factor(Genus)),
# alpha=0.3, size=0.5,
# show.legend = F, width = 0.3, height = 0.1)+
stat_summary(aes(x=Species, y=Affected, fill=Genus),
shape=21,
fun.data = "mean_cl_boot", colour = "black")+
geom_abline(slope = 0, intercept = 0.5, colour="grey")+ facet_grid(~Year_F)
Figure_Count_spp2 
#ggsave(file="Outputs/Figure_2c.svg", plot=Figure_Count_spp2, width=8, height=6.5)
# By sensitivity
Figure_Count_spp2<-ggplot(data=Binomial.data[
Binomial.data$Year!=2014, ]) +
theme(axis.title.x = element_blank(),
legend.title = element_blank())+
#theme(axis.text.x = element_text(angle = 90, vjust = 0.5))+
scale_x_discrete(name="Year")+
scale_y_continuous(name="Proportion of colonies affected",
limits = c(0,1.1),
breaks = c(0, 0.25, 0.5, .75, 1))+
# geom_jitter(aes(x=Year_F, y=Affected,
# colour=as.factor(Sensitivity)),
# alpha=0.3, size=0.5,
# show.legend = F, width = 0.3, height = 0.1)+
stat_summary(aes(x=Year_F, y=Affected, fill=Sensitivity),
shape=21,
fun.data = "mean_cl_boot", colour = "black")
Figure_Count_spp2
Plots by DHW
Spp - Figure S3
Summary_dhw_year<-Binomial.data %>%
group_by(Year) %>%
dplyr::summarise(meanDHW = mean(DHW_max),
seDHW = std.error(DHW_max,na.rm),
sdDHW = sd(DHW_max),
minDHW = min(DHW_max),
maxDHW = max(DHW_max))
Summary_plot_spp<-Binomial.data %>%
group_by(Year, Species) %>%
dplyr::summarise(meanCon = mean(Affected),
sdCon = sd(Affected),
seCon = std.error(Affected,na.rm),
minCon = min(Affected),
maxCon = max(Affected),
n())
Summary_plot_spp<-merge(x = Summary_plot_spp,
y= Summary_dhw_year,
by = c("Year"),
all.x = TRUE)
Summary_plot_sppFigure_spp_dhw<-ggplot(Binomial.data,
aes(x=DHW_max, y=Affected)) +
theme(axis.text.x = element_text(angle = 0, vjust = 0.5))+
geom_abline(slope = 0, intercept = 0.5, colour="grey")+
geom_jitter(data=Binomial.data,
aes(x=DHW_max, y=Affected, colour=Year_F),
width = 0.3, height = 0.05,
alpha=0.3, size=0.5,
show.legend = F)+
stat_summary(aes(x=DHW_max, y=Affected, fill=Year_F),
shape=21, alpha=0.9,
fun.data = "mean_cl_boot", colour = "black")+
scale_x_continuous(name="Max DHW",
#limits = c(0,8),
breaks = c(0, 2, 4, 6, 8))+
scale_y_continuous(name="Proportion of colonies affected",
limits = c(0, 1.1),
breaks = c(0, 0.25, 0.5, .75, 1),
expand = c(0.02, 0.02))+
geom_smooth(span=2, colour="grey", se=F)+
scale_colour_manual(values = c("#FDDBC7", "#67A9CF", "#2166ac", "black",
"#EF8A62", "#B2182B"))+
scale_fill_manual(values = c("#FDDBC7", "#67A9CF", "#2166ac", "black",
"#EF8A62", "#B2182B"))+
facet_wrap(Species~.)
Figure_spp_dhw
#ggsave(file="Outputs/Figure_2014_DHW.svg", plot=Figure_spp_dhw, width=5.5, height=5.5)Figure_spp_year<-ggplot(data=Summary_plot_spp[Summary_plot_spp$Year!=2014, ] ,
aes(x=as.factor(Year), y=meanCon)) +
# scale_x_continuous(name="Max DHW",
# #limits = c(0,8),
# breaks = c(0, 2, 4, 6, 8))+
scale_y_continuous(name="Proportion of colonies affected",
#limits = c(0, 1.1),
breaks = c(0, 0.25, 0.5, .75, 1),
expand = c(0.1, 0.1))+
# binomial_smooth(data=Binomial.data,
# aes(x=DHW_max, y=Affected), colour="black", alpha=0.2)+
geom_jitter(data=Binomial.data[Binomial.data$Year!="2014",],
aes(x=as.factor(Year), y=Affected, colour=as.factor(Year)),
width = 0.3, height = 0.05,
alpha=0.3, size=0.5,
show.legend = F)+
geom_point(shape=21, size=3, aes(fill=as.factor(Year)))+
geom_errorbar(aes(
ymax = meanCon + seCon, ymin = meanCon - seCon),
show.legend = F, width=0.1, colour="black" )+
scale_colour_manual(values = c("#FDDBC7", "#67A9CF", "#2166ac",
"#EF8A62", "#B2182B"))+
scale_fill_manual(values = c("#FDDBC7", "#67A9CF", "#2166ac",
"#EF8A62", "#B2182B"))+
geom_abline(slope = 0, intercept = 0.5, colour="grey")+
facet_wrap(Species~., ncol=3)
Figure_spp_year
#ggsave(file="Outputs/Figure2e.svg", plot=Figure_spp_year, width=5.5, height=5.5)# Figure_spp_dhw_H<-ggplot(data=Binomial.data, aes(x=Affected)) +
# scale_colour_manual(values = c("#FDDBC7", "#67A9CF", "#2166ac",
# "#EF8A62", "#B2182B"))+
# scale_fill_manual(values = c("#FDDBC7", "#67A9CF", "#2166ac",
# "#EF8A62", "#B2182B"))+
# geom_density(aes(fill=Year_F), alpha=0.5, position = "stack")+
# facet_wrap(~Species)
# Figure_spp_dhw_H
#
# Figure_spp_dhw_H<-ggplot(data=Binomial.data, aes(x=Affected)) +
# geom_density(aes(fill=Genus), alpha=0.5, position = "stack")+
# facet_wrap(~Year_F)
# Figure_spp_dhw_H
# Figure_spp_dhw_H<-ggplot(data=Binomial.data, aes(x=Affected)) +
# geom_density(aes(fill=Genus), alpha=0.5, position = "identity")+
# facet_wrap(~Year_F)
# Figure_spp_dhw_HGenera
Summary_plot_genus<-Binomial.data %>%
group_by(Year, Genus) %>%
dplyr::summarise(meanCon = mean(Affected),
sdCon = sd(Affected),
seCon = std.error(Affected,na.rm),
minCon = min(Affected),
maxCon = max(Affected))
Summary_plot_genus<-merge(x = Summary_plot_genus,
y= Summary_dhw_year,
by = c("Year"),
all.x = TRUE)
Summary_plot_genusFigure_gen_dhw<-ggplot(data=Summary_plot_genus,
aes(x=meanDHW, y=meanCon)) +
theme_bw()+
theme(legend.position = "bottom",
panel.grid= element_blank())+
scale_x_continuous(name="Max DHW",
#limits = c(0,8),
breaks = c(0, 2, 4, 6, 8))+
scale_y_continuous(name="Proportion of colonies affected",
limits = c(0, 1),
breaks = c(0, 0.25, 0.5, .75, 1))+
# binomial_smooth(data=Binomial.data,
# aes(x=DHW_max, y=Affected), colour="black", alpha=0.2)+
geom_smooth(span=2, colour="grey", se=F)+
geom_jitter(data=Binomial.data,
aes(x=DHW_max, y=Affected, colour=Year_F),
width = 0.3, height = 0.05,
alpha=0.3, size=0.5, show.legend = F)+
geom_errorbar(aes(
ymax = meanCon + seCon, ymin = meanCon - seCon, colour=as.factor(Year)),
show.legend = F)+
geom_point(shape=21, size=2, aes(fill=as.factor(Year)))+
scale_colour_manual(values = c("#FDDBC7", "#67A9CF", "#2166ac", "blue",
"#EF8A62", "#B2182B"))+
scale_fill_manual(values = c("#FDDBC7", "#67A9CF", "#2166ac", "blue",
"#EF8A62", "#B2182B"))
Figure_gen_dhw + facet_wrap(Genus~.)
#ggsave(file="Outputs/Figure2c.svg", plot=Figure_2c, width=4, height=4)Sensitivity
Summary_plot_sen<-Binomial.data %>%
group_by(Year, Sensitivity) %>%
dplyr::summarise(meanCon = mean(Affected),
sdCon = sd(Affected),
seCon = std.error(Affected,na.rm),
minCon = min(Affected),
maxCon = max(Affected))
Summary_plot_sen<-merge(x = Summary_plot_sen,
y= Summary_dhw_year, by = c("Year"),
all.x = TRUE)
Summary_plot_senFigure_sen_dhw<-ggplot(data=Summary_plot_sen,
aes(x=meanDHW, y=meanCon)) +
scale_x_continuous(name="Max DHW",
#limits = c(0,8),
breaks = c(0, 2, 4, 6, 8))+
scale_y_continuous(name="Proportion of colonies affected",
#limits = c(-.01, 1.1),
breaks = c(0, 0.25, 0.5, .75, 1))+
# binomial_smooth(data=Binomial.data,
# aes(x=DHW_max, y=Affected), colour="black", alpha=0.2)+
geom_smooth(span=2, colour="grey", se=F)+
geom_jitter(data=Binomial.data,
aes(x=DHW_max, y=Affected, colour=Year_F),
width = 0.3, height = 0.05,
alpha=0.3, size=0.5, show.legend = F)+
geom_errorbar(aes(
ymax = meanCon + seCon, ymin = meanCon - seCon))+
geom_point(shape=21, size=2, aes(fill=as.factor(Year)))+
scale_colour_manual(values = c("#FDDBC7", "#67A9CF", "#2166ac", "black",
"#EF8A62", "#B2182B"))+
scale_fill_manual(values = c("#FDDBC7", "#67A9CF", "#2166ac", "black",
"#EF8A62", "#B2182B"))+
facet_wrap(Sensitivity~.)
Figure_sen_dhw
Figure_Count_spp2<-ggplot(data=Summary_plot_sen[Summary_plot_sen$Year!="2014",],
aes(x=as.factor(Year), y=meanCon)) +
theme(axis.title.x = element_blank(),
legend.title = element_blank()
)+
#theme(axis.text.x = element_text(angle = 90, vjust = 0.5))+
scale_x_discrete(name="Year")+
scale_y_continuous(name="Proportion of colonies affected",
limits = c(0,1.1),
breaks = c(0, 0.25, 0.5, .75, 1))+
geom_errorbar(aes(
ymax = meanCon + seCon, ymin = meanCon - seCon),
width = 0.1)+
geom_point(shape=21, size=3, aes(fill=as.factor(Sensitivity)))+
scale_colour_manual(values =
c("#D55E00", "#E69F00", "#009E73"))+
scale_fill_manual(values =
c("#D55E00", "#E69F00", "#009E73"))
Figure_Count_spp2
#ggsave(file="Outputs/Figure_2.svg", plot=Figure_Count_spp2, width=4, height=4)Binomial models
Binomial.data$Sensitivity<-factor(Binomial.data$Sensitivity,
levels = c("Pocillopora",
"Other","Millepora"))# Remove 2014 data with no heat stress
Binomial.data_spp<-Binomial.data[Binomial.data$Year!="2014", ]All_spp_m2<-glm(Affected ~ Year_F * Species,
family = binomial(link = "probit"),
data = Binomial.data_spp)
summary(All_spp_m2)$coefficients## Estimate Std. Error
## (Intercept) 5.575875e+00 143.5344873
## Year_F1997 -5.368799e-09 167.5241612
## Year_F1998 1.108011e+00 0.6279038
## Year_F2015 -7.787395e-09 226.2146746
## Year_F2016 1.313853e-08 580.0877164
## SpeciesPavona varians -3.670548e+00 143.5346487
## SpeciesPavona clavus -6.087159e-09 196.1715039
## SpeciesPocillopora elegans -4.374248e+00 143.5345645
## SpeciesPocillopora damicornis -5.028509e+00 143.5345067
## SpeciesGardineroseris planulata -5.694890e+00 143.5345452
## SpeciesPorites lobata -6.643446e+00 143.5356841
## Year_F1997:SpeciesPavona varians -2.977896e+00 167.5243562
## Year_F1998:SpeciesPavona varians -2.041306e+00 0.6752736
## Year_F2015:SpeciesPavona varians -6.634601e-01 226.2149988
## Year_F2016:SpeciesPavona varians -9.880057e-01 580.0878038
## Year_F1997:SpeciesPavona clavus -6.667495e+00 214.3489798
## Year_F1998:SpeciesPavona clavus -1.108011e+00 201.1113175
## Year_F2015:SpeciesPavona clavus -4.011149e+00 262.7818625
## Year_F2016:SpeciesPavona clavus -4.811165e+00 595.3007146
## Year_F1997:SpeciesPocillopora elegans -2.641158e+00 167.5244857
## Year_F1998:SpeciesPocillopora elegans -3.810724e+00 0.7064904
## Year_F2015:SpeciesPocillopora elegans -2.006223e+00 226.2148010
## Year_F2016:SpeciesPocillopora elegans -1.842294e+00 580.0878040
## Year_F1997:SpeciesPocillopora damicornis -1.445842e+00 167.5242392
## Year_F1998:SpeciesPocillopora damicornis -2.722948e+00 0.6760832
## Year_F2015:SpeciesPocillopora damicornis -8.480755e-01 226.2147521
## Year_F2016:SpeciesPocillopora damicornis -9.028569e-01 580.0877606
## Year_F1997:SpeciesGardineroseris planulata -5.456860e+00 202.1660262
## Year_F1998:SpeciesGardineroseris planulata -1.119014e+00 0.6515428
## Year_F2015:SpeciesGardineroseris planulata 9.951583e-01 226.2149310
## Year_F2016:SpeciesGardineroseris planulata 1.148972e+00 580.0877917
## Year_F1997:SpeciesPorites lobata 5.368841e-09 167.5262120
## Year_F2015:SpeciesPorites lobata 6.972063e-01 226.2155151
## Year_F2016:SpeciesPorites lobata 1.067571e+00 580.0880690
## z value Pr(>|z|)
## (Intercept) 3.884694e-02 9.690124e-01
## Year_F1997 -3.204791e-11 1.000000e+00
## Year_F1998 1.764619e+00 7.762780e-02
## Year_F2015 -3.442480e-11 1.000000e+00
## Year_F2016 2.264921e-11 1.000000e+00
## SpeciesPavona varians -2.557256e-02 9.795983e-01
## SpeciesPavona clavus -3.102978e-11 1.000000e+00
## SpeciesPocillopora elegans -3.047523e-02 9.756881e-01
## SpeciesPocillopora damicornis -3.503345e-02 9.720531e-01
## SpeciesGardineroseris planulata -3.967610e-02 9.683514e-01
## SpeciesPorites lobata -4.628428e-02 9.630837e-01
## Year_F1997:SpeciesPavona varians -1.777590e-02 9.858176e-01
## Year_F1998:SpeciesPavona varians -3.022931e+00 2.503392e-03
## Year_F2015:SpeciesPavona varians -2.932874e-03 9.976599e-01
## Year_F2016:SpeciesPavona varians -1.703200e-03 9.986410e-01
## Year_F1997:SpeciesPavona clavus -3.110579e-02 9.751852e-01
## Year_F1998:SpeciesPavona clavus -5.509441e-03 9.956041e-01
## Year_F2015:SpeciesPavona clavus -1.526418e-02 9.878214e-01
## Year_F2016:SpeciesPavona clavus -8.081908e-03 9.935516e-01
## Year_F1997:SpeciesPocillopora elegans -1.576580e-02 9.874212e-01
## Year_F1998:SpeciesPocillopora elegans -5.393879e+00 6.895274e-08
## Year_F2015:SpeciesPocillopora elegans -8.868664e-03 9.929239e-01
## Year_F2016:SpeciesPocillopora elegans -3.175887e-03 9.974660e-01
## Year_F1997:SpeciesPocillopora damicornis -8.630644e-03 9.931138e-01
## Year_F1998:SpeciesPocillopora damicornis -4.027534e+00 5.636492e-05
## Year_F2015:SpeciesPocillopora damicornis -3.748984e-03 9.970088e-01
## Year_F2016:SpeciesPocillopora damicornis -1.556414e-03 9.987582e-01
## Year_F1997:SpeciesGardineroseris planulata -2.699197e-02 9.784661e-01
## Year_F1998:SpeciesGardineroseris planulata -1.717484e+00 8.589082e-02
## Year_F2015:SpeciesGardineroseris planulata 4.399171e-03 9.964900e-01
## Year_F2016:SpeciesGardineroseris planulata 1.980687e-03 9.984196e-01
## Year_F1997:SpeciesPorites lobata 3.204777e-11 1.000000e+00
## Year_F2015:SpeciesPorites lobata 3.082045e-03 9.975409e-01
## Year_F2016:SpeciesPorites lobata 1.840359e-03 9.985316e-01anova(All_spp_m2)summary(All_spp_m2)##
## Call:
## glm(formula = Affected ~ Year_F * Species, family = binomial(link = "probit"),
## data = Binomial.data_spp)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.66926 -0.55525 0.00016 0.60157 2.32725
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error z value
## (Intercept) 5.576e+00 1.435e+02 0.039
## Year_F1997 -5.369e-09 1.675e+02 0.000
## Year_F1998 1.108e+00 6.279e-01 1.765
## Year_F2015 -7.787e-09 2.262e+02 0.000
## Year_F2016 1.314e-08 5.801e+02 0.000
## SpeciesPavona varians -3.671e+00 1.435e+02 -0.026
## SpeciesPavona clavus -6.087e-09 1.962e+02 0.000
## SpeciesPocillopora elegans -4.374e+00 1.435e+02 -0.030
## SpeciesPocillopora damicornis -5.029e+00 1.435e+02 -0.035
## SpeciesGardineroseris planulata -5.695e+00 1.435e+02 -0.040
## SpeciesPorites lobata -6.643e+00 1.435e+02 -0.046
## Year_F1997:SpeciesPavona varians -2.978e+00 1.675e+02 -0.018
## Year_F1998:SpeciesPavona varians -2.041e+00 6.753e-01 -3.023
## Year_F2015:SpeciesPavona varians -6.635e-01 2.262e+02 -0.003
## Year_F2016:SpeciesPavona varians -9.880e-01 5.801e+02 -0.002
## Year_F1997:SpeciesPavona clavus -6.667e+00 2.143e+02 -0.031
## Year_F1998:SpeciesPavona clavus -1.108e+00 2.011e+02 -0.006
## Year_F2015:SpeciesPavona clavus -4.011e+00 2.628e+02 -0.015
## Year_F2016:SpeciesPavona clavus -4.811e+00 5.953e+02 -0.008
## Year_F1997:SpeciesPocillopora elegans -2.641e+00 1.675e+02 -0.016
## Year_F1998:SpeciesPocillopora elegans -3.811e+00 7.065e-01 -5.394
## Year_F2015:SpeciesPocillopora elegans -2.006e+00 2.262e+02 -0.009
## Year_F2016:SpeciesPocillopora elegans -1.842e+00 5.801e+02 -0.003
## Year_F1997:SpeciesPocillopora damicornis -1.446e+00 1.675e+02 -0.009
## Year_F1998:SpeciesPocillopora damicornis -2.723e+00 6.761e-01 -4.028
## Year_F2015:SpeciesPocillopora damicornis -8.481e-01 2.262e+02 -0.004
## Year_F2016:SpeciesPocillopora damicornis -9.029e-01 5.801e+02 -0.002
## Year_F1997:SpeciesGardineroseris planulata -5.457e+00 2.022e+02 -0.027
## Year_F1998:SpeciesGardineroseris planulata -1.119e+00 6.515e-01 -1.717
## Year_F2015:SpeciesGardineroseris planulata 9.952e-01 2.262e+02 0.004
## Year_F2016:SpeciesGardineroseris planulata 1.149e+00 5.801e+02 0.002
## Year_F1997:SpeciesPorites lobata 5.369e-09 1.675e+02 0.000
## Year_F1998:SpeciesPorites lobata NA NA NA
## Year_F2015:SpeciesPorites lobata 6.972e-01 2.262e+02 0.003
## Year_F2016:SpeciesPorites lobata 1.068e+00 5.801e+02 0.002
## Pr(>|z|)
## (Intercept) 0.9690
## Year_F1997 1.0000
## Year_F1998 0.0776 .
## Year_F2015 1.0000
## Year_F2016 1.0000
## SpeciesPavona varians 0.9796
## SpeciesPavona clavus 1.0000
## SpeciesPocillopora elegans 0.9757
## SpeciesPocillopora damicornis 0.9721
## SpeciesGardineroseris planulata 0.9684
## SpeciesPorites lobata 0.9631
## Year_F1997:SpeciesPavona varians 0.9858
## Year_F1998:SpeciesPavona varians 0.0025 **
## Year_F2015:SpeciesPavona varians 0.9977
## Year_F2016:SpeciesPavona varians 0.9986
## Year_F1997:SpeciesPavona clavus 0.9752
## Year_F1998:SpeciesPavona clavus 0.9956
## Year_F2015:SpeciesPavona clavus 0.9878
## Year_F2016:SpeciesPavona clavus 0.9936
## Year_F1997:SpeciesPocillopora elegans 0.9874
## Year_F1998:SpeciesPocillopora elegans 6.90e-08 ***
## Year_F2015:SpeciesPocillopora elegans 0.9929
## Year_F2016:SpeciesPocillopora elegans 0.9975
## Year_F1997:SpeciesPocillopora damicornis 0.9931
## Year_F1998:SpeciesPocillopora damicornis 5.64e-05 ***
## Year_F2015:SpeciesPocillopora damicornis 0.9970
## Year_F2016:SpeciesPocillopora damicornis 0.9988
## Year_F1997:SpeciesGardineroseris planulata 0.9785
## Year_F1998:SpeciesGardineroseris planulata 0.0859 .
## Year_F2015:SpeciesGardineroseris planulata 0.9965
## Year_F2016:SpeciesGardineroseris planulata 0.9984
## Year_F1997:SpeciesPorites lobata 1.0000
## Year_F1998:SpeciesPorites lobata NA
## Year_F2015:SpeciesPorites lobata 0.9975
## Year_F2016:SpeciesPorites lobata 0.9985
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 3002.3 on 2213 degrees of freedom
## Residual deviance: 1768.1 on 2180 degrees of freedom
## AIC: 1836.1
##
## Number of Fisher Scoring iterations: 16step(All_spp_m2)## Start: AIC=1836.09
## Affected ~ Year_F * Species
##
## Df Deviance AIC
## <none> 1768.1 1836.1
## - Year_F:Species 23 2007.8 2029.8##
## Call: glm(formula = Affected ~ Year_F * Species, family = binomial(link = "probit"),
## data = Binomial.data_spp)
##
## Coefficients:
## (Intercept)
## 5.576e+00
## Year_F1997
## -5.369e-09
## Year_F1998
## 1.108e+00
## Year_F2015
## -7.787e-09
## Year_F2016
## 1.314e-08
## SpeciesPavona varians
## -3.671e+00
## SpeciesPavona clavus
## -6.087e-09
## SpeciesPocillopora elegans
## -4.374e+00
## SpeciesPocillopora damicornis
## -5.029e+00
## SpeciesGardineroseris planulata
## -5.695e+00
## SpeciesPorites lobata
## -6.643e+00
## Year_F1997:SpeciesPavona varians
## -2.978e+00
## Year_F1998:SpeciesPavona varians
## -2.041e+00
## Year_F2015:SpeciesPavona varians
## -6.635e-01
## Year_F2016:SpeciesPavona varians
## -9.880e-01
## Year_F1997:SpeciesPavona clavus
## -6.667e+00
## Year_F1998:SpeciesPavona clavus
## -1.108e+00
## Year_F2015:SpeciesPavona clavus
## -4.011e+00
## Year_F2016:SpeciesPavona clavus
## -4.811e+00
## Year_F1997:SpeciesPocillopora elegans
## -2.641e+00
## Year_F1998:SpeciesPocillopora elegans
## -3.811e+00
## Year_F2015:SpeciesPocillopora elegans
## -2.006e+00
## Year_F2016:SpeciesPocillopora elegans
## -1.842e+00
## Year_F1997:SpeciesPocillopora damicornis
## -1.446e+00
## Year_F1998:SpeciesPocillopora damicornis
## -2.723e+00
## Year_F2015:SpeciesPocillopora damicornis
## -8.481e-01
## Year_F2016:SpeciesPocillopora damicornis
## -9.029e-01
## Year_F1997:SpeciesGardineroseris planulata
## -5.457e+00
## Year_F1998:SpeciesGardineroseris planulata
## -1.119e+00
## Year_F2015:SpeciesGardineroseris planulata
## 9.952e-01
## Year_F2016:SpeciesGardineroseris planulata
## 1.149e+00
## Year_F1997:SpeciesPorites lobata
## 5.369e-09
## Year_F1998:SpeciesPorites lobata
## NA
## Year_F2015:SpeciesPorites lobata
## 6.972e-01
## Year_F2016:SpeciesPorites lobata
## 1.068e+00
##
## Degrees of Freedom: 2213 Total (i.e. Null); 2180 Residual
## Null Deviance: 3002
## Residual Deviance: 1768 AIC: 1836drop1(All_spp_m2, test="LRT")## Check for over/underdispersion in the model
E2 <- resid(All_spp_m2, type = "pearson")
N <- nrow(Binomial.data_spp)
p <- length(coef(All_spp_m2))
sum(E2^2) / (N - p)## [1] 0.8435062All_emmc<-emmeans(All_spp_m2, ~Year_F, by="Species", trans = "log", type = "response")
All_emmc<-multcomp::cld(All_emmc)
All_emmc<-All_emmc[order(All_emmc$Species, All_emmc$Year_F),]
All_emmcwrite.csv(All_emmc, "Outputs/Species_comparissions.csv", row.names = F)
All_emmc_1<-emmeans(All_spp_m2, ~Species, by="Year_F", trans = "log", type = "response")
All_emmc_1<-multcomp::cld(All_emmc_1)
All_emmc_1<-All_emmc_1[order(All_emmc_1$Year_F, All_emmc_1$Species),]
All_emmc_1write.csv(All_emmc_1, "Outputs/year_comparissions.csv", row.names = F)2. Supplementary data - other locations
Data
Number of colonies sampled
Summary_plot_Non_uva<-Binomial_other %>%
group_by(Year, Species) %>%
dplyr::summarise(meanCon = mean(Affected),
sdCon = sd(Affected),
seCon = std.error(Affected,na.rm),
minCon = min(Affected),
maxCon = max(Affected),
n())
Summary_plot_Non_uvaPlots by Year
Facet species - Figure S4
binomial_smooth <- function(...) {
geom_smooth(method = "glm", method.args = list(family = "binomial"), ...)
}
Figure_Count<-ggplot(data=Binomial_other,
aes(x=Year_F, y=Affected)) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5))+
scale_x_discrete(name="Year")+
scale_y_continuous(name="Proportion of colonies affected",
limits = c(0,1.1),
breaks = c(0, 0.25, 0.5, .75, 1))+
geom_abline(slope = 0, intercept = 0.5, colour="grey")+
geom_jitter(aes(x=Year_F, y=Affected,
colour=as.factor(Location)),
show.legend = F,
alpha=0.5, size=0.8, width = 0.2, height = 0.1)+
stat_summary(aes(x=Year_F, y=Affected), alpha=0.9,
fun.data = "mean_cl_boot", colour = "black")
Figure_Count<-Figure_Count + facet_wrap(~Species, ncol=3)
Figure_Count
#ggsave(file="Outputs/Figure2_otherLocations.svg", plot=Figure_Count, width=5, height=5.5)Facet years
Figure_Count_spp2_o<-ggplot(data=Binomial_other) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
axis.title.x = element_blank(),
legend.title = element_blank())+
scale_x_discrete(name="Year")+
scale_y_continuous(name="Proportion of colonies affected",
limits = c(0,1.1),
breaks = c(0, 0.25, 0.5, .75, 1))+
geom_jitter(aes(x=Species, y=Affected,
colour=as.factor(Genus)),
alpha=0.3, size=0.5,
show.legend = F, width = 0.3, height = 0.1)+
stat_summary(aes(x=Species, y=Affected, fill=Genus),
shape=21,
fun.data = "mean_cl_boot", colour = "black")+
geom_abline(slope = 0, intercept = 0.5, colour="grey")+
facet_grid(Location~Year_F)
Figure_Count_spp2_o 
#ggsave(file="Outputs/Figure_2c.svg", plot=Figure_Count_spp2, width=8, height=6.5)
Figure_Count_spp2_o<-ggplot(data=Binomial_other) +
theme(axis.title.x = element_blank(),
legend.title = element_blank())+
#theme(axis.text.x = element_text(angle = 90, vjust = 0.5))+
scale_x_discrete(name="Year")+
scale_y_continuous(name="Proportion of colonies affected",
limits = c(0,1.1),
breaks = c(0, 0.25, 0.5, .75, 1))+
# geom_jitter(aes(x=Year_F, y=Affected,
# colour=as.factor(Sensitivity)),
# alpha=0.3, size=0.5,
# show.legend = F, width = 0.3, height = 0.1)+
stat_summary(aes(x=Year_F, y=Affected, fill=Sensitivity),
shape=21,
fun.data = "mean_cl_boot", colour = "black")
Figure_Count_spp2_o
Auguie, Baptiste. 2017. GridExtra: Miscellaneous Functions for "Grid" Graphics. https://CRAN.R-project.org/package=gridExtra.
Genz, Alan, and Frank Bretz. 2009. Computation of Multivariate Normal and T Probabilities. Lecture Notes in Statistics. Heidelberg: Springer-Verlag.
Genz, Alan, Frank Bretz, Tetsuhisa Miwa, Xuefei Mi, and Torsten Hothorn. 2021. Mvtnorm: Multivariate Normal and T Distributions. http://mvtnorm.R-forge.R-project.org.
Hothorn, Torsten. 2019. TH.data: TH’s Data Archive. https://CRAN.R-project.org/package=TH.data.
Hothorn, Torsten, Frank Bretz, and Peter Westfall. 2008. “Simultaneous Inference in General Parametric Models.” Biometrical Journal 50 (3): 346–63.
———. 2021. Multcomp: Simultaneous Inference in General Parametric Models. https://CRAN.R-project.org/package=multcomp.
Lenth, Russell V. 2022. Emmeans: Estimated Marginal Means, Aka Least-Squares Means. https://github.com/rvlenth/emmeans.
Pinheiro, José, Douglas Bates, and R-core. 2021. Nlme: Linear and Nonlinear Mixed Effects Models. https://svn.r-project.org/R-packages/trunk/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/.
Ripley, Brian. 2021. MASS: Support Functions and Datasets for Venables and Ripley’s Mass. http://www.stats.ox.ac.uk/pub/MASS4/.
Terry M. Therneau, and Patricia M. Grambsch. 2000. Modeling Survival Data: Extending the Cox Model. New York: Springer.
Therneau, Terry M. 2021. Survival: Survival Analysis. https://github.com/therneau/survival.
Venables, W. N., and B. D. Ripley. 2002. Modern Applied Statistics with S. Fourth. New York: Springer. https://www.stats.ox.ac.uk/pub/MASS4/.
Wickham, Hadley. 2016. Ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York. https://ggplot2.tidyverse.org.
Wickham, Hadley, Winston Chang, Lionel Henry, Thomas Lin Pedersen, Kohske Takahashi, Claus Wilke, Kara Woo, Hiroaki Yutani, and Dewey Dunnington. 2021. 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. 2021. Dplyr: A Grammar of Data Manipulation. https://CRAN.R-project.org/package=dplyr.
Wickham, Hadley, and Dana Seidel. 2020. Scales: Scale Functions for Visualization. https://CRAN.R-project.org/package=scales.