Uva Island bleaching 2
Ana M. Palacio-Castro
2022-07-29
Setup
Import data for binary approach
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_2016.csv",
header = T)
# nrow(Bi_2014_16)
# summary(Bi_2014_16)
Binomial.data<-rbind(Bi_1983_97, Bi_2014_16)
summary(Binomial.data)## Location Date Species
## Uva :2561 03/20/83:1248 Pocillopora damicornis :517
## Canal de Afuera: 52 03/13/98: 604 Porites panamensis :428
## Secas : 234 08/19/15: 119 Pavona varians :392
## 08/16/15: 99 Pocillopora elegans :309
## 04/20/16: 92 Gardineroseris planulata:285
## 08/18/15: 90 Psammocora stellata :213
## (Other) : 595 (Other) :703
## Condition Site
## Affected:1754 Length:2847
## Healthy :1093 Class :character
## Mode :character
##
##
##
## #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.data$Date<-as.Date(Binomial.data$Date,
format="%m/%d/%y")
Binomial.data$Year<-year(Binomial.data$Date)
Binomial.data$Year_F<-as.factor(Binomial.data$Year)
Binomial.data<-merge(Binomial.data, Species, by="Species", all.x=TRUE)
Binomial.data$Sensitivity<-factor(Binomial.data$Sensitivity,
levels = c("Millepora",
"Other", "Pocillopora"))
Binomial.data$Genus<-factor(Binomial.data$Genus,
levels = c("Millepora",
"Pavona", "Gardineroseris",
"Porites",
"Psammocora", "Pocillopora"))
Binomial.data$Species<-factor(Binomial.data$Species, levels = c(
"Millepora platyphylla", "Millepora intricata",
"Pavona chiriquiensis", "Pavona varians", "Pavona clavus",
"Porites panamensis", "Porites lobata", "Pavona gigantea", "Psammocora stellata" , "Gardineroseris planulata", "Pocillopora damicornis", "Pocillopora elegans", "Pocillopora"))
Binomial.data$Affected[Binomial.data$Condition=="Healthy"]<-0
Binomial.data$Affected[Binomial.data$Condition=="Affected"]<-1
# filter data
# Remove 2018 - Recovery
Binomial.data<-Binomial.data[Binomial.data$Year!=2018,]
# Remove - species with a few years
Binomial.data<-Binomial.data[Binomial.data$Species!="Millepora platyphylla",]
Binomial.data<-Binomial.data[Binomial.data$Species!="Pavona chiriquiensis",]
Binomial.data<-Binomial.data[Binomial.data$Species!="Psammocora stellata",]
# Remove - non-uva locations
Binomial.data<-Binomial.data[Binomial.data$Location=="Uva",]
Binomial.data<-(droplevels(Binomial.data))
summary(Binomial.data)## Species Location Date Condition
## Pocillopora damicornis :454 Uva:2278 Min. :1983-03-20 Affected:1571
## Porites panamensis :419 1st Qu.:1983-03-20 Healthy : 707
## Pavona varians :368 Median :1998-03-13
## Gardineroseris planulata:277 Mean :1996-01-14
## Pocillopora elegans :251 3rd Qu.:2014-08-24
## Pavona clavus :169 Max. :2016-04-21
## (Other) :340
## Site Year Year_F Genus
## Length:2278 Min. :1983 1983:1086 Millepora : 96
## Class :character 1st Qu.:1983 1998: 543 Pavona :660
## Mode :character Median :1998 2014: 103 Gardineroseris:277
## Mean :1996 2015: 335 Porites :540
## 3rd Qu.:2014 2016: 211 Pocillopora :705
## Max. :2016
##
## Sensitivity Affected
## Millepora : 96 Min. :0.0000
## Other :1477 1st Qu.:0.0000
## Pocillopora: 705 Median :1.0000
## Mean :0.6896
## 3rd Qu.:1.0000
## Max. :1.0000
## Heat stress
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(-(30: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': 2278 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/ 9 levels "Millepora intricata",..: 7 7 2 7 7 1 2 2 7 7 ...
## $ Condition : Factor w/ 2 levels "Affected","Healthy": 1 2 1 2 2 1 1 1 1 2 ...
## $ Site : chr NA NA NA NA ...
## $ Year : num 1983 1983 1983 1983 1983 ...
## $ Year_F : Factor w/ 5 levels "1983","1998",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Genus : Factor w/ 5 levels "Millepora","Pavona",..: 3 3 2 3 3 1 2 2 3 3 ...
## $ Sensitivity: Factor w/ 3 levels "Millepora","Other",..: 2 2 2 2 2 1 2 2 2 2 ...
## $ Affected : num 1 0 1 0 0 1 1 1 1 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 ...Plots by Year
binomial_smooth <- function(...) {
geom_smooth(method = "glm", method.args = list(family = "binomial"), ...)
}
Figure_Count<-ggplot(data=Binomial.data,
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_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,
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 + facet_wrap(Species~.)
Figure_Count + facet_wrap(Genus~.)
Figure_Count + facet_wrap(Sensitivity~.)
#ggsave(file="Outputs/Figure2d.svg", plot=Figure_2d, width=5, height=5)Figure_Count_spp2<-ggplot(data=Binomial.data) +
#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
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(data=Summary_plot_spp,
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(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),
expand = c(0.1, 0.1))+
# binomial_smooth(data=Binomial.data,
# aes(x=DHW_max, y=Affected), colour="black", alpha=0.2)+
geom_smooth(span=2, 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)+
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",
"#EF8A62", "#B2182B"))+
scale_fill_manual(values = c("#FDDBC7", "#67A9CF", "#2166ac",
"#EF8A62", "#B2182B"))
Figure_spp_dhw + facet_wrap(Species~.)
#ggsave(file="Outputs/Figure2c.svg", plot=Figure_2c, width=4, height=4)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_H
Genera
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(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")+
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",
"#EF8A62", "#B2182B"))+
scale_fill_manual(values = c("#FDDBC7", "#67A9CF", "#2166ac",
"#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(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",
"#EF8A62", "#B2182B"))+
scale_fill_manual(values = c("#FDDBC7", "#67A9CF", "#2166ac",
"#EF8A62", "#B2182B"))+
facet_wrap(Sensitivity~.)
Figure_sen_dhw
#ggsave(file="Outputs/Figure2c.svg", plot=Figure_2c, width=4, height=4)Binomial models
Binomial.data$Sensitivity<-factor(Binomial.data$Sensitivity,
levels = c("Pocillopora",
"Other","Millepora"))Thought time
library(jtools)By sensitivity - Pocillopora
Pocillopora.data<- Binomial.data[Binomial.data$Sensitivity=="Pocillopora",]
Pocillopora_m0<-glm(Affected ~ DHW_max,
family="binomial",
data = Pocillopora.data[Pocillopora.data$Year!="2014",])
Pocillopora_m1<-glm(Affected ~ DHW_max*Year,
family="binomial",
data = Pocillopora.data[Pocillopora.data$Year!="2014",])
Pocillopora_m2<-glm(Affected ~ Year,
family="binomial",
data = Pocillopora.data[Pocillopora.data$Year!="2014",])
anova(Pocillopora_m0, Pocillopora_m1, Pocillopora_m2,
test = "Chisq")step(Pocillopora_m1)## Start: AIC=760.39
## Affected ~ DHW_max * Year
##
## Df Deviance AIC
## - DHW_max:Year 1 752.91 758.91
## <none> 752.39 760.39
##
## Step: AIC=758.91
## Affected ~ DHW_max + Year
##
## Df Deviance AIC
## <none> 752.91 758.91
## - DHW_max 1 819.28 823.28
## - Year 1 911.16 915.16##
## Call: glm(formula = Affected ~ DHW_max + Year, family = "binomial",
## data = Pocillopora.data[Pocillopora.data$Year != "2014",
## ])
##
## Coefficients:
## (Intercept) DHW_max Year
## 147.04725 0.61534 -0.07547
##
## Degrees of Freedom: 694 Total (i.e. Null); 692 Residual
## Null Deviance: 952
## Residual Deviance: 752.9 AIC: 758.9New.data <- expand.grid(
Year = seq(from = min(Binomial.data$Year),
to = max(Binomial.data$Year), by = 2),
DHW_max = seq(0, 9, by=0.5))
fits <- predict(Pocillopora_m1, newdata = New.data, se.fit = TRUE, type="response")
# Get odds from modrl
New.data$prediction <- exp(fits$fit)
New.data$upper <- exp(fits$fit + 1.96 * fits$se.fit)
New.data$lower <- exp(fits$fit - 1.96 * fits$se.fit)
# Convert odds to probabilities
New.data$prediction <- New.data$prediction / (1 + New.data$prediction)
New.data$upper <- New.data$upper / (1 + New.data$upper)
New.data$lower <- New.data$lower / (1 + New.data$lower)
# Plot probabilities
ggplot(New.data, aes(DHW_max, prediction)) +
#geom_errorbar(aes(ymin = lower, ymax = upper, colour = Year),
# width = 0.25, size = 0.5, position = position_dodge(width = 0.4))+
geom_point(aes(fill = Year), shape = 21, size = 3,
position = position_dodge(width = 0.4)) +
geom_line(aes(DHW_max, prediction, group = Year))+
#facet_grid(~education) +
scale_y_continuous(name = "Probability of bleaching/Mortality",
limits = c(0, 1),
labels = scales::percent)
plot_model(Pocillopora_m1, type = "pred",
terms = c("DHW_max", "Year"))
By sensitivity - Millepora
Millepora.data<- Binomial.data[Binomial.data$Sensitivity=="Millepora",]
Millepora_m0<-glm(Affected ~ DHW_max,
family="binomial",
data = Millepora.data)
Millepora_m1<-glm(Affected ~ DHW_max * Year,
family="binomial",
data = Millepora.data)
Millepora_m2<-glm(Affected ~ Year,
family="binomial",
data = Millepora.data)
anova(Millepora_m0, Millepora_m1, Millepora_m2, test = "Chisq")step(Millepora_m1)## Start: AIC=15.48
## Affected ~ DHW_max * Year
##
## Df Deviance AIC
## - DHW_max:Year 1 7.4813 13.481
## <none> 7.4813 15.481
##
## Step: AIC=13.48
## Affected ~ DHW_max + Year
##
## Df Deviance AIC
## - Year 1 7.481 11.481
## <none> 7.481 13.481
## - DHW_max 1 64.642 68.642
##
## Step: AIC=11.48
## Affected ~ DHW_max
##
## Df Deviance AIC
## <none> 7.481 11.481
## - DHW_max 1 83.213 85.213##
## Call: glm(formula = Affected ~ DHW_max, family = "binomial", data = Millepora.data)
##
## Coefficients:
## (Intercept) DHW_max
## -2.708 4.050
##
## Degrees of Freedom: 95 Total (i.e. Null); 94 Residual
## Null Deviance: 83.21
## Residual Deviance: 7.481 AIC: 11.48anova(Millepora_m1)New.data <- expand.grid(
Year = seq(from = min(Millepora.data$Year),
to = max(Millepora.data$Year), by = 2),
DHW_max = seq(0, 9, by=0.5))
fits <- predict(Millepora_m0, newdata = New.data, se.fit = TRUE, type="response")
# Get odds from modrl
New.data$prediction <- exp(fits$fit)
New.data$upper <- exp(fits$fit + 1.96 * fits$se.fit)
New.data$lower <- exp(fits$fit - 1.96 * fits$se.fit)
# Convert odds to probabilities
New.data$prediction <- New.data$prediction / (1 + New.data$prediction)
New.data$upper <- New.data$upper / (1 + New.data$upper)
New.data$lower <- New.data$lower / (1 + New.data$lower)
# Plot probabilities
ggplot(New.data, aes(DHW_max, prediction)) +
#geom_errorbar(aes(ymin = lower, ymax = upper, colour = Year),
# width = 0.25, size = 0.5, position = position_dodge(width = 0.4))+
geom_point(aes(fill = Year), shape = 21, size = 3,
position = position_dodge(width = 0.4)) +
geom_line(aes(DHW_max, prediction, colour = Year))+
#facet_grid(~education) +
scale_y_continuous(name = "Probability of bleaching/Mortality",
limits = c(0, 1),
labels = scales::percent)
plot_model(Millepora_m0, type = "pred",
terms = c("DHW_max"))
By sensitivity - Others
Other.data<- Binomial.data[Binomial.data$Sensitivity=="Other",]
Other_m0<-glm(Affected ~ DHW_max,
family="binomial",
data = Other.data)
Other_m1<-glm(Affected ~ DHW_max*Year,
family="binomial",
data = Other.data)
Other_m2<-glm(Affected ~ Year,
family="binomial",
data = Other.data)
anova(Other_m0, Other_m1, Other_m1, test = "Chisq")step(Other_m1)## Start: AIC=1477.26
## Affected ~ DHW_max * Year
##
## Df Deviance AIC
## <none> 1469.3 1477.3
## - DHW_max:Year 1 1509.6 1515.6##
## Call: glm(formula = Affected ~ DHW_max * Year, family = "binomial",
## data = Other.data)
##
## Coefficients:
## (Intercept) DHW_max Year DHW_max:Year
## 394.91239 -49.40740 -0.19724 0.02476
##
## Degrees of Freedom: 1476 Total (i.e. Null); 1473 Residual
## Null Deviance: 1682
## Residual Deviance: 1469 AIC: 1477anova(Other_m1, test = "Chisq" )New.data <- expand.grid(
Year = seq(from = min(Other.data$Year),
to = max(Other.data$Year), by = 2),
DHW_max = seq(0, 9, by=0.5))
fits <- predict(Other_m1, newdata = New.data, se.fit = TRUE, type="response")
# Get odds from modrl
New.data$prediction <- exp(fits$fit)
New.data$upper <- exp(fits$fit + 1.96 * fits$se.fit)
New.data$lower <- exp(fits$fit - 1.96 * fits$se.fit)
# Convert odds to probabilities
New.data$prediction <- New.data$prediction / (1 + New.data$prediction)
New.data$upper <- New.data$upper / (1 + New.data$upper)
New.data$lower <- New.data$lower / (1 + New.data$lower)
# Plot probabilities
ggplot(New.data, aes(DHW_max, prediction)) +
#geom_errorbar(aes(ymin = lower, ymax = upper, colour = Year),
# width = 0.25, size = 0.5, position = position_dodge(width = 0.4))+
geom_point(aes(fill = Year), shape = 21, size = 3,
position = position_dodge(width = 0.4)) +
geom_line(aes(DHW_max, prediction, colour = Year))+
#facet_grid(~education) +
scale_y_continuous(name = "Probability of bleaching/Mortality",
limits = c(0, 1),
labels = scales::percent)
plot_model(Other_m1, type = "pred",
terms = c("DHW_max", "Year"))
All
Model_m0<-glm(Affected ~ DHW_max * Sensitivity,
family="binomial",
data = Binomial.data)
Model_m1<-glm(Affected ~ DHW_max + Year * Sensitivity,
family="binomial",
data =Binomial.data)
Model_m2<-glm(Affected ~ Year_F * Sensitivity,
family="binomial",
data = Binomial.data)
anova(Model_m0, Model_m1, Model_m2, test = "Chisq")step(Model_m1)## Start: AIC=2329.85
## Affected ~ DHW_max + Year * Sensitivity
##
## Df Deviance AIC
## <none> 2315.8 2329.8
## - Year:Sensitivity 2 2324.0 2334.0
## - DHW_max 1 2466.9 2478.9##
## Call: glm(formula = Affected ~ DHW_max + Year * Sensitivity, family = "binomial",
## data = Binomial.data)
##
## Coefficients:
## (Intercept) DHW_max
## 140.26844 0.30311
## Year SensitivityOther
## -0.07114 -40.22061
## SensitivityMillepora Year:SensitivityOther
## 112.15610 0.02092
## Year:SensitivityMillepora
## -0.05435
##
## Degrees of Freedom: 2277 Total (i.e. Null); 2271 Residual
## Null Deviance: 2822
## Residual Deviance: 2316 AIC: 2330anova(Model_m1)New.data <- expand.grid(
Sensitivity = unique(Binomial.data$Sensitivity),
Year = seq(from = min(Binomial.data$Year),
to = max(Binomial.data$Year), by = 2),
DHW_max = seq(0, 9, by=0.5))
fits <- predict(Model_m1, newdata = New.data,
se.fit = TRUE, type="response")
# Get odds from modrl
New.data$prediction <- exp(fits$fit)
New.data$upper <- exp(fits$fit + 1.96 * fits$se.fit)
New.data$lower <- exp(fits$fit - 1.96 * fits$se.fit)
# Convert odds to probabilities
New.data$prediction <- New.data$prediction / (1 + New.data$prediction)
New.data$upper <- New.data$upper / (1 + New.data$upper)
New.data$lower <- New.data$lower / (1 + New.data$lower)
# Plot probabilities
ggplot(New.data, aes(DHW_max, prediction)) +
#geom_errorbar(aes(ymin = lower, ymax = upper, colour = Year),
# width = 0.25, size = 0.5, position = position_dodge(width = 0.4))+
geom_point(aes(fill = Year), shape = 21, size = 3,
position = position_dodge(width = 0.4)) +
geom_line(aes(DHW_max, prediction, colour = Year))+
facet_wrap(~Sensitivity) +
scale_y_continuous(name = "Probability of bleaching/Mortality",
limits = c(0, 1),
labels = scales::percent)
plot_model(Model_m1, type = "pred",
terms = c("Year", "Sensitivity"))
plot_model(Model_m1, type = "pred",
terms = c("DHW_max", "Sensitivity"))
plot_model(Model_m1, type = "pred",
terms = c("DHW_max", "Sensitivity", "Year"))
All species ENSO
ENSOs.data<- Binomial.data[Binomial.data$Year!="2014",]
ENSOs_m0<-glm(Affected ~ Year_F,
family="binomial",
data = ENSOs.data)
ENSOs_m1<-glm(Affected ~ Sensitivity,
family="binomial",
data = ENSOs.data)
ENSOs_m2<-glm(Affected ~ Year_F * Sensitivity,
family="binomial",
data = ENSOs.data)
ENSOs_m3<-glm(Affected ~ DHW_max + Year_F * Sensitivity,
family="binomial",
data = ENSOs.data)
ENSOs_m4<-glm(Affected ~ DHW_max * Sensitivity,
family="binomial",
data = ENSOs.data)
anova(ENSOs_m0, ENSOs_m1, ENSOs_m2, test = "Chisq")anova(ENSOs_m2, ENSOs_m3, test = "Chisq")anova(ENSOs_m2, ENSOs_m4, test = "Chisq")New.data <- expand.grid(
Sensitivity = unique(ENSOs.data$Sensitivity),
# DHW_max = seq(0, 9, by=0.5),
Year_F = unique(ENSOs.data$Year_F))
fits <- predict(ENSOs_m2, newdata = New.data,
se.fit = TRUE, type="response")
# Get odds from modrl
New.data$prediction <- exp(fits$fit)
New.data$upper <- exp(fits$fit + 1.96 * fits$se.fit)
New.data$lower <- exp(fits$fit - 1.96 * fits$se.fit)
# Convert odds to probabilities
New.data$prediction <- New.data$prediction / (1 + New.data$prediction)
New.data$upper <- New.data$upper / (1 + New.data$upper)
New.data$lower <- New.data$lower / (1 + New.data$lower)
# Plot probabilities (1)
ggplot(New.data, aes(Sensitivity, prediction)) +
geom_errorbar(aes(ymin = lower,
ymax = upper,
colour = Year_F),
width = 0.25, size = 0.5,
position = position_dodge(width = 0.4))+
geom_point(aes(fill = Year_F), shape = 21, size = 3,
position = position_dodge(width = 0.4)) +
#facet_grid(~education) +
scale_y_continuous(name = "Probability of bleaching/Mortality",
limits = c(0, 1),
labels = scales::percent)
# Plot probabilities (2)
ggplot(New.data, aes(Year_F, prediction)) +
geom_errorbar(aes(ymin = lower,
ymax = upper,
colour = Sensitivity),
width = 0.25, size = 0.5,
position = position_dodge(width = 0.4))+
geom_point(aes(fill = Sensitivity), shape = 21, size = 3,
position = position_dodge(width = 0.4)) +
#facet_grid(~education) +
scale_y_continuous(name = "Probability of bleaching/Mortality",
limits = c(0, 1),
labels = scales::percent)
POst-hoc
anova(ENSOs_m2, test = "Chisq")summary(ENSOs_m2)##
## Call:
## glm(formula = Affected ~ Year_F * Sensitivity, family = "binomial",
## data = ENSOs.data)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.0162 -0.4673 0.5300 0.7967 2.1301
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.1387 0.1116 10.203 < 2e-16 ***
## Year_F1998 -3.2982 0.3693 -8.931 < 2e-16 ***
## Year_F2015 -2.0116 0.2354 -8.545 < 2e-16 ***
## Year_F2016 -1.8831 0.3001 -6.274 3.52e-10 ***
## SensitivityOther 0.7534 0.1644 4.583 4.58e-06 ***
## SensitivityMillepora 15.4274 353.7936 0.044 0.96522
## Year_F1998:SensitivityOther 2.3910 0.4025 5.940 2.85e-09 ***
## Year_F2015:SensitivityOther 0.9324 0.3073 3.034 0.00241 **
## Year_F2016:SensitivityOther 0.7659 0.3684 2.079 0.03763 *
## Year_F1998:SensitivityMillepora NA NA NA NA
## Year_F2015:SensitivityMillepora 2.0116 557.5894 0.004 0.99712
## Year_F2016:SensitivityMillepora 1.8831 1429.8397 0.001 0.99895
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2585.0 on 2174 degrees of freedom
## Residual deviance: 2176.7 on 2164 degrees of freedom
## AIC: 2198.7
##
## Number of Fisher Scoring iterations: 15effect_plot(ENSOs_m2,
pred = Year_F, interval = TRUE,
#partial.residuals = TRUE,
jitter = c(0.5, 0),
y.label = "% proportion affected")
effect_plot(ENSOs_m2,
pred = Sensitivity, interval = TRUE,
#partial.residuals = TRUE,
jitter = c(0.5, 0),
y.label = "% proportion affected")
ENSOs_m3<-glm(Affected ~ DHW_max + Year_F * Sensitivity,
family="binomial",
data = ENSOs.data)
New.data <- expand.grid(
DHW_max = seq(3, 8.5, by=0.5),
Sensitivity = unique(ENSOs.data$Sensitivity),
Year_F= unique(ENSOs.data$Year_F))
fits <- predict(ENSOs_m3, newdata = New.data,
se.fit = TRUE, type="response")
# Get odds from modrl
New.data$prediction <- exp(fits$fit)
New.data$upper <- exp(fits$fit + 1.96 * fits$se.fit)
New.data$lower <- exp(fits$fit - 1.96 * fits$se.fit)
# Convert odds to probabilities
New.data$prediction <- New.data$prediction / (1 + New.data$prediction)
New.data$upper <- New.data$upper / (1 + New.data$upper)
New.data$lower <- New.data$lower / (1 + New.data$lower)
# Plot probabilities (1)
ggplot(New.data, aes(DHW_max, prediction)) +
geom_errorbar(aes(ymin = lower,
ymax = upper,
colour = Year_F),
width = 0.25, size = 0.5,
position = position_dodge(width = 0.4))+
geom_point(aes(fill = Year_F), shape = 21, size = 3,
position = position_dodge(width = 0.4)) +
facet_grid(~Sensitivity) +
scale_y_continuous(name = "Probability of bleaching/Mortality",
limits = c(0, 1),
labels = scales::percent)
# Plot probabilities (2)
ggplot(New.data, aes(DHW_max, prediction)) +
geom_errorbar(aes(ymin = lower,
ymax = upper,
colour = Sensitivity),
width = 0.25, size = 0.5,
position = position_dodge(width = 0.4))+
geom_point(aes(fill = Sensitivity), shape = 21, size = 3,
position = position_dodge(width = 0.4)) +
facet_grid(~Year_F) +
scale_y_continuous(name = "Probability of bleaching/Mortality",
limits = c(0, 1),
labels = scales::percent)
effect_plot(ENSOs_m3,
pred = Year_F, interval = TRUE,
#partial.residuals = TRUE,
jitter = c(0.5, 0),
y.label = "% proportion affected")
effect_plot(ENSOs_m3,
pred = Sensitivity, interval = TRUE,
#partial.residuals = TRUE,
jitter = c(0.5, 0),
y.label = "% proportion affected")
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. 2020. 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. 2021. 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. 2020. 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.