PET bleaching
Ana M. Palacio-Castro
2022-04-06
Setup
Coral health data
Import data
# Read data
Bleaching_proportion <- read.csv(
"Data/CoralHealth-ETP-Master_1808-4Ana.csv", header=TRUE, sep = ",")
summary(Bleaching_proportion)## Year Mo Date Time
## Min. :2012 Min. : 3.00 04/20/16: 192 : 457
## 1st Qu.:2015 1st Qu.: 4.00 08/19/15: 173 08/22/15 09:00 AM: 124
## Median :2015 Median : 8.00 08/16/15: 157 02:00:00 PM : 122
## Mean :2015 Mean : 6.39 03/12/17: 140 08/19/15 10:00 AM: 119
## 3rd Qu.:2016 3rd Qu.: 8.00 08/22/15: 135 08/16/15 10:00 AM: 113
## Max. :2018 Max. :10.00 04/21/16: 133 04/21/16 09:30 AM: 98
## (Other) :1870 (Other) :1767
## Region Island Site Location
## Eastern Panama : 241 Uva :1074 Reef :918 Ridley : 529
## Northern Galapagos: 508 Secas-Dog : 379 Uva-Ridley:335 Forereef : 479
## Southern Galapagos: 205 Darwin : 294 Uva-Reef :284 Millepora: 348
## Western Panama :1846 Uva-Ridley: 245 Dog :192 Dog : 162
## Wolf : 214 Uva-reef :105 Reef : 160
## Floreana : 205 Mainland :104 4x5 : 120
## (Other) : 389 (Other) :862 (Other) :1002
## Transect Z.Category Z..m. Z..ft.
## UvRf-0-6-1: 123 0-6 :1007 Min. : 0.30 : 674
## ICA : 98 15-24: 702 1st Qu.: 4.30 66 : 122
## UvRf-0-6-2: 94 25 : 179 Median :10.00 14 : 92
## SaRf-0-6-1: 91 7-14 : 912 Mean :10.78 36 : 83
## T1 : 82 3rd Qu.:15.80 8 : 83
## T4 : 71 Max. :40.00 15 : 75
## (Other) :2241 NA's :2 (Other):1671
## Tr..Type Tr..Length Observer Morph
## Ind. Colonies :1977 :2292 Ana Palacio: 12 :2546
## Permanent Transect : 378 i : 189 JM : 147 br : 45
## Haphazard Intercept: 176 10 : 172 PF : 481 En : 1
## ICA : 83 10x1 : 36 RB : 278 encr: 35
## Permanent Intercept: 79 4 : 21 TBS :1178 ms : 167
## Transect : 36 30 : 19 VB : 704 nd : 5
## (Other) : 71 (Other): 71 plt : 1
## Spp Spp2 Length Width Height
## POLO :387 POCSP :627 Min. : 1.00 30 : 130 Min. : 0.00
## PAGI :384 POLO :387 1st Qu.: 15.00 10 : 121 1st Qu.: 8.00
## MIIN :363 PAGI :384 Median : 28.00 15 : 116 Median : 14.00
## PODA :360 MIIN :363 Mean : 43.99 12 : 105 Mean : 23.42
## PACL :353 PACL :353 3rd Qu.: 52.25 20 : 102 3rd Qu.: 30.00
## PAVA :246 PAVA :246 Max. :650.00 14 : 91 Max. :300.00
## (Other):707 (Other):440 NA's :4 (Other):2135 NA's :9
## X..Bleach X..Pale Old Recent
## Min. : 1.0 Min. : 1.00 Min. : 1.00 Min. : 1.00
## 1st Qu.: 8.0 1st Qu.: 10.00 1st Qu.:10.00 1st Qu.: 3.00
## Median : 20.0 Median : 30.00 Median :20.00 Median : 5.00
## Mean : 32.9 Mean : 34.22 Mean :26.38 Mean : 13.55
## 3rd Qu.: 50.0 3rd Qu.: 50.00 3rd Qu.:40.00 3rd Qu.: 10.00
## Max. :100.0 Max. :100.00 Max. :98.00 Max. :100.00
## NA's :2059 NA's :1678 NA's :768 NA's :2287
## Bites Pred X.Damsel Damsel.spp Lithoph
## :2364 Min. : 1.000 :2660 :2665 :2615
## LOTS : 62 1st Qu.: 3.000 1 : 119 Sbeeb : 42 LOTS : 42
## 10 : 53 Median : 5.000 2 : 8 Sacap : 38 4 : 20
## 5 : 37 Mean : 5.923 3 : 7 Stbe : 25 1 : 18
## 4 : 34 3rd Qu.: 7.750 5 : 2 SACAP : 13 5 : 18
## 3 : 31 Max. :20.000 7 : 1 Sflav : 4 6 : 18
## (Other): 219 NA's :2774 PF: 3 (Other): 13 (Other): 69
## Pink X.Algae Alg.spp X.FLCY
## :2627 :1854 :1857 :2608
## 4 : 33 2 : 174 FIL : 620 2 : 45
## 2 : 31 4 : 154 PEYS : 91 4 : 40
## 5 : 26 6 : 106 FLCY : 71 6 : 19
## 3 : 21 5 : 89 Peyssonnellia: 25 10 : 14
## 6 : 14 10 : 88 CAUL : 24 8 : 12
## (Other): 48 (Other): 335 (Other) : 112 (Other): 62
## X..Barnacles X..Worm Eucidaris X..Snail X..Disease
## :2615 :2728 Min. : 1.000 :2786 Min. : 5
## LOTS : 45 15 : 10 1st Qu.: 1.000 1 : 5 1st Qu.:10
## 2 : 25 10 : 9 Median : 2.000 2 : 3 Median :10
## 4 : 22 2 : 7 Mean : 2.562 4 : 2 Mean :14
## 6 : 14 20 : 7 3rd Qu.: 3.000 0 : 1 3rd Qu.:10
## 5 : 13 1 : 6 Max. :22.000 10 : 1 Max. :35
## (Other): 66 (Other): 33 NA's :2736 (Other): 2 NA's :2795
## Disease.Type Black.Sponge Sponge COTS
## Mode:logical Min. : 3.000 Min. : 1 Min. :1
## NA's:2800 1st Qu.: 5.000 1st Qu.: 1 1st Qu.:1
## Median : 8.000 Median : 1 Median :1
## Mean : 7.455 Mean : 4 Mean :1
## 3rd Qu.:10.000 3rd Qu.: 2 3rd Qu.:1
## Max. :15.000 Max. :15 Max. :1
## NA's :2789 NA's :2795 NA's :2799
## Notes BL.0
## :2094 Min. : 0.000
## many bites : 21 1st Qu.: 0.000
## No full date on sheet, but month and year : 21 Median : 0.000
## Not totally sure if this was ICA or transect: 14 Mean : 8.707
## SED 1 : 14 3rd Qu.: 2.000
## 4SED : 13 Max. :100.000
## (Other) : 623
## PA.0 Old.0 Recent.0 X
## Min. : 0.00 Min. : 0.00 Min. : 0.000 Min. :0.0000
## 1st Qu.: 0.00 1st Qu.: 0.00 1st Qu.: 0.000 1st Qu.:0.0000
## Median : 0.00 Median :10.00 Median : 0.000 Median :0.0000
## Mean : 13.71 Mean :19.15 Mean : 2.482 Mean :0.2646
## 3rd Qu.: 20.00 3rd Qu.:30.00 3rd Qu.: 0.000 3rd Qu.:1.0000
## Max. :100.00 Max. :98.00 Max. :100.000 Max. :1.0000
## #Bleaching_proportion<-Bleaching_proportion %>% drop_na(value)
Bleaching_proportion$Date<-as.Date(Bleaching_proportion$Date,
format="%m/%d/%y")
#Bleaching_proportion$Year<-year(Bleaching_proportion$Date)
Bleaching_proportion$Year_F<-as.factor(Bleaching_proportion$Year)
Bleaching_proportion$Z.Category<-factor(Bleaching_proportion$Z.Category, levels =
c("0-6", "7-14", "15-24", "25"))
Bleaching_WP<-Bleaching_proportion[Bleaching_proportion$Region=="Western Panama", ]
# Remove spp
Bleaching_WP_SPP<-Bleaching_WP[Bleaching_WP$Spp2!="PACH", ] # Pavona chiriquiensis
Bleaching_WP_SPP<-Bleaching_WP_SPP[Bleaching_WP_SPP$Spp2!="PSST", ] #Psammocora stellata
Bleaching_WP_SPP<-Bleaching_WP_SPP[Bleaching_WP_SPP$Spp2!="PSSP", ] #Psammocora species
Bleaching_WP_SPP<-Bleaching_WP_SPP[Bleaching_WP_SPP$Spp2!="PSSU", ] #Psammocora superficialis
Bleaching_WP_SPP<-droplevels(Bleaching_WP_SPP)
Bleaching_WP_SPP$Spp2<-factor(Bleaching_WP_SPP$Spp2, levels = c("MIIN", "GAPL", "PACL", "PAGI", "PAVA",
"POLO", "POPA", "POCSP" ))
Bleaching_WP_SPP$Spp<-factor(Bleaching_WP_SPP$Spp, levels = c("PODA", "MIIN", "POSP", "GAPL", "PACL",
"POEL", "PAGI", "PAVA",
"POEY", "POPA", "POLO"))
# Remove 2014 and 2018
Bleaching_WP_SPP_2015_16<-Bleaching_WP_SPP[Bleaching_WP_SPP$Year!="2014", ] # 2014
Bleaching_WP_SPP_2015_16<-Bleaching_WP_SPP_2015_16[Bleaching_WP_SPP_2015_16$Year!="2018", ] #2018
summary(Bleaching_WP_SPP_2015_16$Spp2)## MIIN GAPL PACL PAGI PAVA POLO POPA POCSP
## 242 72 120 131 137 127 116 454 summary(Bleaching_WP_SPP_2015_16$Spp)## PODA MIIN POSP GAPL PACL POEL PAGI PAVA POEY POPA POLO
## 265 242 0 72 120 183 131 137 6 116 127 #Bleaching_WP_SPP<-Bleaching_WP_SPP %>% drop_na(Z.Category)Bleaching
All locations
Plot_all_sites <- ggplot(Bleaching_proportion, aes(
x=as.factor(Year_F), y=BL.0)) +
geom_boxplot(outlier.shape = NA) +
stat_boxplot(geom = 'errorbar')+
geom_jitter(aes(fill=Spp2), shape=21, alpha=0.8, width = 0.15) +
#stat_summary(fun.data = mean_se, geom = "errorbar")+
#stat_summary(fun=mean, geom="point", size=2) +
#stat_summary(fun="mean_cl_normal", geom = "errorbar", size=2) +
#stat_summary(fun.y=median, geom="point", size=2, shape=14) +
facet_grid(~Region) + theme(legend.position = "right") +
xlab("Year") + ylab("Bleaching percentaje")
Plot_all_sites
#ggsave(file="Outputs/Bleaching2.svg", plot=Plot_all_boxplot2, width=6, height=5)Species and years (Only western Panama)
Plot_all_years <- ggplot(Bleaching_WP_SPP, aes(
x=as.factor(Year_F), y=BL.0)) +
geom_boxplot(outlier.shape = NA) +
stat_boxplot(geom = 'errorbar')+
geom_jitter(aes(fill=Spp), shape=21, alpha=0.8, width = 0.15) +
#stat_summary(fun.data = mean_se, geom = "errorbar")+
#stat_summary(fun=mean, geom="point", size=2) +
#stat_summary(fun="mean_cl_normal", geom = "errorbar", size=2) +
#stat_summary(fun.y=median, geom="point", size=2, shape=14) +
facet_grid(Spp2~Z.Category) + theme(legend.position = "right") +
xlab("Year") + ylab("Bleaching percentaje")
Plot_all_years
#ggsave(file="Outputs/Bleaching2.svg", plot=Plot_all_boxplot2, width=6, height=5)Species and years (Only western Panama, 2015 and 2016)
Option 1
Plot_2015_16 <- ggplot(Bleaching_WP_SPP_2015_16, aes(
x=as.factor(Year_F), y=BL.0)) +
geom_boxplot(outlier.shape = NA) +
stat_boxplot(geom = 'errorbar')+
geom_jitter(aes(fill=Spp), shape=21, alpha=0.8, width = 0.20) +
#stat_summary(fun.data = mean_se, geom = "errorbar")+
#stat_summary(fun=mean, geom="point", size=2) +
#stat_summary(fun="mean_cl_normal", geom = "errorbar", size=2) +
#stat_summary(fun.y=median, geom="point", size=2, shape=14) +
facet_grid(Spp2~Z.Category) + theme(legend.position = "right") +
xlab("Year") + ylab("Bleaching percentaje")
Plot_2015_16
#ggsave(file="Outputs/Bleaching2.svg", plot=Plot_all_boxplot2, width=6, height=5)Option 2
Plot_2015_16_b <- ggplot(Bleaching_WP_SPP_2015_16, aes(
x=as.factor(Z.Category), y=BL.0)) +
geom_boxplot(outlier.shape = NA) +
stat_boxplot(geom = 'errorbar')+
geom_jitter(aes(fill=Spp), shape=21, alpha=0.8, width = 0.20) +
#stat_summary(fun.data = mean_se, geom = "errorbar")+
#stat_summary(fun=mean, geom="point", size=2) +
#stat_summary(fun="mean_cl_normal", geom = "errorbar", size=2) +
#stat_summary(fun.y=median, geom="point", size=2, shape=14) +
facet_grid(Spp2~Year_F) + theme(legend.position = "right") +
xlab("Year") + ylab("Bleaching percentaje")
Plot_2015_16_b
#ggsave(file="Outputs/Bleaching2.svg", plot=Plot_all_boxplot2, width=6, height=5)Option 3
binomial_smooth <- function(...) {
geom_smooth(method = "glm", method.args = list(family = "binomial"), ...)
}
Bleaching_b <- ggplot(Bleaching_WP_SPP, aes(
x=Z..m., y=BL.0/100)) +
binomial_smooth()+
geom_jitter(aes(fill=Spp), shape=21, alpha=0.8, width = 0.20) +
#stat_summary(fun.data = mean_se, geom = "errorbar")+
#stat_summary(fun=mean, geom="point", size=2) +
#stat_summary(fun="mean_cl_normal", geom = "errorbar", size=2) +
#stat_summary(fun.y=median, geom="point", size=2, shape=14) +
facet_grid(Spp2~Year_F) + theme(legend.position = "right") +
xlab("Depth (m)") + ylab("Bleaching proportion")
Bleaching_b
Option 4
Bleaching_binomial <- ggplot(Bleaching_WP_SPP_2015_16, aes(
x=Z..m., y=BL.0/100), colour=Spp2) +
binomial_smooth(aes(colour=Spp2), se=F)+
geom_jitter(aes(fill=Spp2), shape=21, alpha=0.8, width = 0.20) +
#stat_summary(fun.data = mean_se, geom = "errorbar")+
#stat_summary(fun=mean, geom="point", size=2) +
#stat_summary(fun="mean_cl_normal", geom = "errorbar", size=2) +
#stat_summary(fun.y=median, geom="point", size=2, shape=14) +
facet_grid(~Year_F) + theme(legend.position = "right") +
xlab("Depth (m)") + ylab("Bleaching proportion")
Bleaching_binomial
Paling
All locations
Pale_all_sites <- ggplot(Bleaching_proportion, aes(
x=as.factor(Year_F), y=PA.0)) +
geom_boxplot(outlier.shape = NA) +
stat_boxplot(geom = 'errorbar')+
geom_jitter(aes(fill=Spp2), shape=21, alpha=0.8, width = 0.15) +
#stat_summary(fun.data = mean_se, geom = "errorbar")+
#stat_summary(fun=mean, geom="point", size=2) +
#stat_summary(fun="mean_cl_normal", geom = "errorbar", size=2) +
#stat_summary(fun.y=median, geom="point", size=2, shape=14) +
facet_grid(~Region) + theme(legend.position = "right") +
xlab("Year") + ylab("Paling percentaje")
Pale_all_sites
#ggsave(file="Outputs/Bleaching2.svg", plot=Plot_all_boxplot2, width=6, height=5)Species and years (Only western Panama)
Plot_all_years <- ggplot(Bleaching_WP_SPP, aes(
x=as.factor(Year_F), y=PA.0)) +
geom_boxplot(outlier.shape = NA) +
stat_boxplot(geom = 'errorbar')+
geom_jitter(aes(fill=Spp), shape=21, alpha=0.8, width = 0.15) +
#stat_summary(fun.data = mean_se, geom = "errorbar")+
#stat_summary(fun=mean, geom="point", size=2) +
#stat_summary(fun="mean_cl_normal", geom = "errorbar", size=2) +
#stat_summary(fun.y=median, geom="point", size=2, shape=14) +
facet_grid(Spp2~Z.Category) + theme(legend.position = "right") +
xlab("Year") + ylab("Paling percentaje")
Plot_all_years
#ggsave(file="Outputs/Bleaching2.svg", plot=Plot_all_boxplot2, width=6, height=5)Species and years (Only western Panama, 2015 and 2016)
Option 1
Pale_2015_16 <- ggplot(Bleaching_WP_SPP_2015_16, aes(
x=as.factor(Year_F), y=PA.0)) +
geom_boxplot(outlier.shape = NA) +
stat_boxplot(geom = 'errorbar')+
geom_jitter(aes(fill=Spp), shape=21, alpha=0.8, width = 0.20) +
#stat_summary(fun.data = mean_se, geom = "errorbar")+
#stat_summary(fun=mean, geom="point", size=2) +
#stat_summary(fun="mean_cl_normal", geom = "errorbar", size=2) +
#stat_summary(fun.y=median, geom="point", size=2, shape=14) +
facet_grid(Spp2~Z.Category) + theme(legend.position = "right") +
xlab("Year") + ylab("Paling percentaje")
Pale_2015_16
Option 2
Pale_2015_16b <- ggplot(Bleaching_WP_SPP_2015_16, aes(
x=as.factor(Z.Category), y=PA.0)) +
geom_boxplot(outlier.shape = NA) +
stat_boxplot(geom = 'errorbar')+
geom_jitter(aes(fill=Spp), shape=21, alpha=0.8, width = 0.20) +
#stat_summary(fun.data = mean_se, geom = "errorbar")+
#stat_summary(fun=mean, geom="point", size=2) +
#stat_summary(fun="mean_cl_normal", geom = "errorbar", size=2) +
#stat_summary(fun.y=median, geom="point", size=2, shape=14) +
facet_grid(Spp2~Year_F) + theme(legend.position = "right") +
xlab("Year") + ylab("Paling percentaje")
Pale_2015_16b
#ggsave(file="Outputs/Bleaching2.svg", plot=Plot_all_boxplot2, width=6, height=5)Option 3
Pale_2015_16_c <- ggplot(Bleaching_WP_SPP, aes(
x=Z..m., y=PA.0/100)) +
binomial_smooth()+
geom_jitter(aes(fill=Spp), shape=21, alpha=0.8, width = 0.20) +
#stat_summary(fun.data = mean_se, geom = "errorbar")+
#stat_summary(fun=mean, geom="point", size=2) +
#stat_summary(fun="mean_cl_normal", geom = "errorbar", size=2) +
#stat_summary(fun.y=median, geom="point", size=2, shape=14) +
facet_grid(Spp2~Year_F) + theme(legend.position = "right") +
xlab("Depth (m)") + ylab("Paling proportion")
Pale_2015_16_c
Option 4
Pale_binomial <- ggplot(Bleaching_WP_SPP_2015_16, aes(
x=Z..m., y=PA.0/100), colour=Spp2) +
binomial_smooth(aes(colour=Spp2), se=F)+
geom_jitter(aes(fill=Spp2), shape=21, alpha=0.8, width = 0.20) +
#stat_summary(fun.data = mean_se, geom = "errorbar")+
#stat_summary(fun=mean, geom="point", size=2) +
#stat_summary(fun="mean_cl_normal", geom = "errorbar", size=2) +
#stat_summary(fun.y=median, geom="point", size=2, shape=14) +
facet_grid(~Year_F) + theme(legend.position = "right") +
xlab("Depth (m)") + ylab("Bleaching proportion")
Pale_binomial
Paling + Bleaching
All locations
Plot_all_sites <- ggplot(Bleaching_proportion, aes(
x=as.factor(Year_F), y=BL.0+PA.0)) +
geom_boxplot(outlier.shape = NA) +
stat_boxplot(geom = 'errorbar')+
geom_jitter(aes(fill=Spp2), shape=21, alpha=0.8, width = 0.15) +
#stat_summary(fun.data = mean_se, geom = "errorbar")+
#stat_summary(fun=mean, geom="point", size=2) +
#stat_summary(fun="mean_cl_normal", geom = "errorbar", size=2) +
#stat_summary(fun.y=median, geom="point", size=2, shape=14) +
facet_grid(~Region) + theme(legend.position = "right") +
xlab("Year") + ylab("Bleaching percentaje")
Plot_all_sites
#ggsave(file="Outputs/Bleaching2.svg", plot=Plot_all_boxplot2, width=6, height=5)Species and years (Only western Panama)
Plot_all_years <- ggplot(Bleaching_WP_SPP, aes(
x=as.factor(Year_F), y=(BL.0+PA.0))) +
geom_boxplot(outlier.shape = NA) +
stat_boxplot(geom = 'errorbar')+
geom_jitter(aes(fill=Spp), shape=21, alpha=0.8, width = 0.15) +
#stat_summary(fun.data = mean_se, geom = "errorbar")+
#stat_summary(fun=mean, geom="point", size=2) +
#stat_summary(fun="mean_cl_normal", geom = "errorbar", size=2) +
#stat_summary(fun.y=median, geom="point", size=2, shape=14) +
facet_grid(Spp2~Z.Category) + theme(legend.position = "right") +
xlab("Year") + ylab("Bleaching + paling ")
Plot_all_years
#ggsave(file="Outputs/Bleaching2.svg", plot=Plot_all_boxplot2, width=6, height=5)Species and years (Only western Panama, 2015 and 2016)
Option 1
Pale_2015_16 <- ggplot(Bleaching_WP_SPP_2015_16, aes(
x=as.factor(Year_F), y=(BL.0+PA.0))) +
geom_boxplot(outlier.shape = NA) +
stat_boxplot(geom = 'errorbar')+
geom_jitter(aes(fill=Spp), shape=21, alpha=0.8, width = 0.20) +
#stat_summary(fun.data = mean_se, geom = "errorbar")+
#stat_summary(fun=mean, geom="point", size=2) +
#stat_summary(fun="mean_cl_normal", geom = "errorbar", size=2) +
#stat_summary(fun.y=median, geom="point", size=2, shape=14) +
facet_grid(Spp2~Z.Category) + theme(legend.position = "right") +
xlab("Year") + ylab("Bleaching + paling ")
Pale_2015_16
#ggsave(file="Outputs/Bleaching2.svg", plot=Plot_all_boxplot2, width=6, height=5)Option 2
Pale_2015_16b <- ggplot(Bleaching_WP_SPP_2015_16, aes(
x=as.factor(Z.Category), y=(BL.0+PA.0))) +
geom_boxplot(outlier.shape = NA) +
stat_boxplot(geom = 'errorbar')+
geom_jitter(aes(fill=Spp), shape=21, alpha=0.8, width = 0.20) +
#stat_summary(fun.data = mean_se, geom = "errorbar")+
#stat_summary(fun=mean, geom="point", size=2) +
#stat_summary(fun="mean_cl_normal", geom = "errorbar", size=2) +
#stat_summary(fun.y=median, geom="point", size=2, shape=14) +
facet_grid(Spp2~Year_F) + theme(legend.position = "right") +
xlab("Year") + ylab("Bleaching + paling ")
Pale_2015_16b
#ggsave(file="Outputs/Bleaching2.svg", plot=Plot_all_boxplot2, width=6, height=5)Option 3
Plot_2015_16_c <- ggplot(Bleaching_WP_SPP, aes(
x=Z..m., y=(BL.0+PA.0)/100)) +
binomial_smooth()+
geom_jitter(aes(fill=Spp), shape=21, alpha=0.8, width = 0.20) +
#stat_summary(fun.data = mean_se, geom = "errorbar")+
#stat_summary(fun=mean, geom="point", size=2) +
#stat_summary(fun="mean_cl_normal", geom = "errorbar", size=2) +
#stat_summary(fun.y=median, geom="point", size=2, shape=14) +
facet_grid(Spp2~Year_F) + theme(legend.position = "right") +
xlab("Depth (m)") + ylab("Bleaching + paling ")
Plot_2015_16_c
Option 4
2 data points are percentajes > 100%
PaleB_binomial <- ggplot(Bleaching_WP_SPP_2015_16, aes(
x=Z..m., y=(PA.0+BL.0)/100), colour=Spp2) +
binomial_smooth(aes(colour=Spp2), se=F)+
geom_jitter(aes(fill=Spp2), shape=21, alpha=0.8, width = 0.20) +
#stat_summary(fun.data = mean_se, geom = "errorbar")+
#stat_summary(fun=mean, geom="point", size=2) +
#stat_summary(fun="mean_cl_normal", geom = "errorbar", size=2) +
#stat_summary(fun.y=median, geom="point", size=2, shape=14) +
facet_grid(~Year_F) + theme(legend.position = "right") +
xlab("Depth (m)") + ylab("Bleaching proportion")
PaleB_binomial
Recent mortality
All locations
Plot_all_recent <- ggplot(Bleaching_proportion, aes(
x=as.factor(Year_F), y=Recent.0)) +
geom_boxplot(outlier.shape = NA) +
stat_boxplot(geom = 'errorbar')+
geom_jitter(aes(fill=Spp2), shape=21, alpha=0.8, width = 0.15) +
#stat_summary(fun.data = mean_se, geom = "errorbar")+
#stat_summary(fun=mean, geom="point", size=2) +
#stat_summary(fun="mean_cl_normal", geom = "errorbar", size=2) +
#stat_summary(fun.y=median, geom="point", size=2, shape=14) +
facet_grid(~Region) + theme(legend.position = "right") +
xlab("Year") + ylab("Bleaching percentaje")
Plot_all_recent
Species and years (Only western Panama)
Plot_recent_years <- ggplot(Bleaching_WP_SPP, aes(
x=as.factor(Year_F), y=(Recent.0))) +
geom_boxplot(outlier.shape = NA) +
stat_boxplot(geom = 'errorbar')+
geom_jitter(aes(fill=Spp), shape=21, alpha=0.8, width = 0.15) +
#stat_summary(fun.data = mean_se, geom = "errorbar")+
#stat_summary(fun=mean, geom="point", size=2) +
#stat_summary(fun="mean_cl_normal", geom = "errorbar", size=2) +
#stat_summary(fun.y=median, geom="point", size=2, shape=14) +
facet_grid(Spp2~Z.Category) + theme(legend.position = "right") +
xlab("Year") + ylab("Recent mortality")
Plot_recent_years
#ggsave(file="Outputs/Bleaching2.svg", plot=Plot_all_boxplot2, width=6, height=5)Species and years (Only western Panama, all years?)
Option 1
Recent_2014_18b <- ggplot(Bleaching_WP_SPP, aes(
x=as.factor(Z.Category), y=(Recent.0))) +
geom_boxplot(outlier.shape = NA) +
stat_boxplot(geom = 'errorbar')+
geom_jitter(aes(fill=Spp), shape=21, alpha=0.8, width = 0.20) +
#stat_summary(fun.data = mean_se, geom = "errorbar")+
#stat_summary(fun=mean, geom="point", size=2) +
#stat_summary(fun="mean_cl_normal", geom = "errorbar", size=2) +
#stat_summary(fun.y=median, geom="point", size=2, shape=14) +
facet_grid(Spp2~Year_F, scales = "free_y") + theme(legend.position = "right") +
xlab("Year") + ylab("Recent mortality")
Recent_2014_18b
Option 2
BS for all spp except M. intrincata 2015 (best data for this particuar spp)
Recent_2014_18c <- ggplot(Bleaching_WP_SPP, aes(
x=Z..m., y=(Recent.0)/100)) +
binomial_smooth()+
geom_jitter(aes(fill=Spp), shape=21, alpha=0.8, width = 0.20) +
#stat_summary(fun.data = mean_se, geom = "errorbar")+
#stat_summary(fun=mean, geom="point", size=2) +
#stat_summary(fun="mean_cl_normal", geom = "errorbar", size=2) +
#stat_summary(fun.y=median, geom="point", size=2, shape=14) +
facet_grid(Spp2~Year_F) + theme(legend.position = "right") +
xlab("Depth (m)") + ylab("Recent mortality")
Recent_2014_18c
Option 4
Recent_binomial <- ggplot(Bleaching_WP_SPP_2015_16, aes(
x=Z..m., y=Recent.0/100), colour=Spp2) +
binomial_smooth(aes(colour=Spp2), se=F)+
geom_jitter(aes(fill=Spp2), shape=21, alpha=0.8, width = 0.20) +
#stat_summary(fun.data = mean_se, geom = "errorbar")+
#stat_summary(fun=mean, geom="point", size=2) +
#stat_summary(fun="mean_cl_normal", geom = "errorbar", size=2) +
#stat_summary(fun.y=median, geom="point", size=2, shape=14) +
facet_grid(~Year_F) + theme(legend.position = "right") +
xlab("Depth (m)") + ylab("Recent mortality")
Recent_binomial
Old mortality
All locations
Plot_all_sites <- ggplot(Bleaching_proportion, aes(
x=as.factor(Year_F), y=Old.0)) +
geom_boxplot(outlier.shape = NA) +
stat_boxplot(geom = 'errorbar')+
geom_jitter(aes(fill=Spp2), shape=21, alpha=0.8, width = 0.15) +
#stat_summary(fun.data = mean_se, geom = "errorbar")+
#stat_summary(fun=mean, geom="point", size=2) +
#stat_summary(fun="mean_cl_normal", geom = "errorbar", size=2) +
#stat_summary(fun.y=median, geom="point", size=2, shape=14) +
facet_grid(~Region) + theme(legend.position = "right") +
xlab("Year") + ylab("Bleaching percentaje")
Plot_all_sites
#ggsave(file="Outputs/Bleaching2.svg", plot=Plot_all_boxplot2, width=6, height=5)Species and years (Only western Panama)
Plot_all_years <- ggplot(Bleaching_WP_SPP, aes(
x=as.factor(Year_F), y=(Old.0))) +
geom_boxplot(outlier.shape = NA) +
stat_boxplot(geom = 'errorbar')+
geom_jitter(aes(fill=Spp), shape=21, alpha=0.8, width = 0.15) +
#stat_summary(fun.data = mean_se, geom = "errorbar")+
#stat_summary(fun=mean, geom="point", size=2) +
#stat_summary(fun="mean_cl_normal", geom = "errorbar", size=2) +
#stat_summary(fun.y=median, geom="point", size=2, shape=14) +
facet_grid(Spp2~Z.Category) + theme(legend.position = "right") +
xlab("Year") + ylab("Old mortality ")
Plot_all_years
#ggsave(file="Outputs/Bleaching2.svg", plot=Plot_all_boxplot2, width=6, height=5)Species and years (Only western Panama, 2015 and 2016)
Option 1
old_2015_16 <- ggplot(Bleaching_WP_SPP_2015_16, aes(
x=as.factor(Year_F), y=(Old.0))) +
geom_boxplot(outlier.shape = NA) +
stat_boxplot(geom = 'errorbar')+
geom_jitter(aes(fill=Spp), shape=21, alpha=0.8, width = 0.20) +
#stat_summary(fun.data = mean_se, geom = "errorbar")+
#stat_summary(fun=mean, geom="point", size=2) +
#stat_summary(fun="mean_cl_normal", geom = "errorbar", size=2) +
#stat_summary(fun.y=median, geom="point", size=2, shape=14) +
facet_grid(Spp2~Z.Category) + theme(legend.position = "right") +
xlab("Year") + ylab("Old mortality ")
old_2015_16
#ggsave(file="Outputs/Bleaching2.svg", plot=Plot_all_boxplot2, width=6, height=5)Option 2
old_2015_16b <- ggplot(Bleaching_WP_SPP_2015_16, aes(
x=as.factor(Z.Category), y=(Old.0))) +
geom_boxplot(outlier.shape = NA) +
stat_boxplot(geom = 'errorbar')+
geom_jitter(aes(fill=Spp), shape=21, alpha=0.8, width = 0.20) +
#stat_summary(fun.data = mean_se, geom = "errorbar")+
#stat_summary(fun=mean, geom="point", size=2) +
#stat_summary(fun="mean_cl_normal", geom = "errorbar", size=2) +
#stat_summary(fun.y=median, geom="point", size=2, shape=14) +
facet_grid(Spp2~Year_F) + theme(legend.position = "right") +
xlab("Year") + ylab("Old mortality ")
old_2015_16b
#ggsave(file="Outputs/Bleaching2.svg", plot=Plot_all_boxplot2, width=6, height=5)Option 3
old_2015_16c <- ggplot(Bleaching_WP_SPP, aes(
x=Z..m., y=(Old.0)/100)) +
binomial_smooth()+
geom_jitter(aes(fill=Spp), shape=21, alpha=0.8, width = 0.20) +
#stat_summary(fun.data = mean_se, geom = "errorbar")+
#stat_summary(fun=mean, geom="point", size=2) +
#stat_summary(fun="mean_cl_normal", geom = "errorbar", size=2) +
#stat_summary(fun.y=median, geom="point", size=2, shape=14) +
facet_grid(Spp2~Year_F) + theme(legend.position = "right") +
xlab("Depht (m)") + ylab("Old mortality ")
old_2015_16c
Models
Bleaching
# Initial Proportion Model by rank
Bleaching_0<-glm (BL.0/100 ~ Z..m., family = binomial,
data=Bleaching_WP_SPP_2015_16, na.action=na.omit)
Bleaching_1<-glm (BL.0/100 ~ Z..m. + Spp2, family = binomial,
data=Bleaching_WP_SPP_2015_16, na.action=na.omit)
Bleaching_2.1<-glm (BL.0/100 ~ Z..m. + Spp2 + Year_F, family = binomial,
data=Bleaching_WP_SPP_2015_16, na.action=na.omit)
Bleaching_2<-glm (BL.0/100 ~ Z..m. * Spp2 + Year_F, family = binomial,
data=Bleaching_WP_SPP_2015_16, na.action=na.omit)
anova(Bleaching_0, Bleaching_1, test = "Chisq") anova(Bleaching_1, Bleaching_2, test = "Chisq") anova(Bleaching_2.1, Bleaching_2, test = "Chisq") step(Bleaching_2)## Start: AIC=839.59
## BL.0/100 ~ Z..m. * Spp2 + Year_F
##
## Df Deviance AIC
## <none> 502.30 839.59
## - Z..m.:Spp2 7 517.24 840.53
## - Year_F 1 549.29 884.59##
## Call: glm(formula = BL.0/100 ~ Z..m. * Spp2 + Year_F, family = binomial,
## data = Bleaching_WP_SPP_2015_16, na.action = na.omit)
##
## Coefficients:
## (Intercept) Z..m. Spp2GAPL Spp2PACL
## -0.196780 -0.129078 1.185667 0.870426
## Spp2PAGI Spp2PAVA Spp2POLO Spp2POPA
## -0.675272 1.100065 0.264444 0.502828
## Spp2POCSP Year_F2016 Z..m.:Spp2GAPL Z..m.:Spp2PACL
## -1.457295 -1.234540 -0.104420 -0.018143
## Z..m.:Spp2PAGI Z..m.:Spp2PAVA Z..m.:Spp2POLO Z..m.:Spp2POPA
## 0.008401 -0.080914 -0.187906 0.032990
## Z..m.:Spp2POCSP
## 0.097479
##
## Degrees of Freedom: 1398 Total (i.e. Null); 1382 Residual
## Null Deviance: 728.7
## Residual Deviance: 502.3 AIC: 839.6 summary(Bleaching_2)##
## Call:
## glm(formula = BL.0/100 ~ Z..m. * Spp2 + Year_F, family = binomial,
## data = Bleaching_WP_SPP_2015_16, na.action = na.omit)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.35359 -0.45899 -0.27006 0.01777 2.54545
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.196780 0.537348 -0.366 0.7142
## Z..m. -0.129078 0.044600 -2.894 0.0038 **
## Spp2GAPL 1.185667 0.824299 1.438 0.1503
## Spp2PACL 0.870426 0.647467 1.344 0.1788
## Spp2PAGI -0.675272 0.685690 -0.985 0.3247
## Spp2PAVA 1.100065 0.643704 1.709 0.0875 .
## Spp2POLO 0.264444 0.746491 0.354 0.7232
## Spp2POPA 0.502828 0.677125 0.743 0.4577
## Spp2POCSP -1.457295 0.628456 -2.319 0.0204 *
## Year_F2016 -1.234540 0.191564 -6.445 1.16e-10 ***
## Z..m.:Spp2GAPL -0.104420 0.118546 -0.881 0.3784
## Z..m.:Spp2PACL -0.018143 0.065739 -0.276 0.7826
## Z..m.:Spp2PAGI 0.008401 0.079031 0.106 0.9153
## Z..m.:Spp2PAVA -0.080914 0.069807 -1.159 0.2464
## Z..m.:Spp2POLO -0.187906 0.127699 -1.471 0.1412
## Z..m.:Spp2POPA 0.032990 0.055194 0.598 0.5500
## Z..m.:Spp2POCSP 0.097479 0.057346 1.700 0.0892 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 728.7 on 1398 degrees of freedom
## Residual deviance: 502.3 on 1382 degrees of freedom
## AIC: 839.59
##
## Number of Fisher Scoring iterations: 6 scatter.smooth(fitted(Bleaching_2), resid(Bleaching_2)); abline(h=0, lty=2)
confint(Bleaching_2) ## CIs using profiled log-likelihood## 2.5 % 97.5 %
## (Intercept) -1.29631335 0.83100613
## Z..m. -0.21762540 -0.04135951
## Spp2GAPL -0.37415340 2.87096542
## Spp2PACL -0.37582888 2.17608962
## Spp2PAGI -2.00699871 0.69956182
## Spp2PAVA -0.13692779 2.40048631
## Spp2POLO -1.17717321 1.76097571
## Spp2POPA -0.80777554 1.86026883
## Spp2POCSP -2.67698501 -0.19731203
## Year_F2016 -1.61950432 -0.86711913
## Z..m.:Spp2GAPL -0.36549027 0.10277142
## Z..m.:Spp2PACL -0.15326299 0.10685254
## Z..m.:Spp2PAGI -0.17640755 0.14785650
## Z..m.:Spp2PAVA -0.23112577 0.04771559
## Z..m.:Spp2POLO -0.46544126 0.03795467
## Z..m.:Spp2POPA -0.07685685 0.14063244
## Z..m.:Spp2POCSP -0.01579201 0.20990151 confint.default(Bleaching_2) ## CIs using standard errors## 2.5 % 97.5 %
## (Intercept) -1.24996299 0.85640207
## Z..m. -0.21649216 -0.04166359
## Spp2GAPL -0.42992989 2.80126347
## Spp2PACL -0.39858521 2.13943689
## Spp2PAGI -2.01919966 0.66865558
## Spp2PAVA -0.16157192 2.36170146
## Spp2POLO -1.19865137 1.72753961
## Spp2POPA -0.82431304 1.82996908
## Spp2POCSP -2.68904552 -0.22554375
## Year_F2016 -1.60999884 -0.85908100
## Z..m.:Spp2GAPL -0.33676576 0.12792563
## Z..m.:Spp2PACL -0.14699009 0.11070313
## Z..m.:Spp2PAGI -0.14649623 0.16329855
## Z..m.:Spp2PAVA -0.21773319 0.05590495
## Z..m.:Spp2POLO -0.43819066 0.06237960
## Z..m.:Spp2POPA -0.07518969 0.14116870
## Z..m.:Spp2POCSP -0.01491778 0.20987537 overdisp_fun <- function(model) {
rdf <- df.residual(model)
rp <- residuals(model,type="deviance")
Pearson.chisq <- sum(rp^2)
prat <- Pearson.chisq/rdf
pval <- pchisq(Pearson.chisq, df=rdf, lower.tail=FALSE)
c(chisq=Pearson.chisq,ratio=prat,rdf=rdf,p=pval)
}
overdisp_fun(Bleaching_2)## chisq ratio rdf p
## 502.2967484 0.3634564 1382.0000000 1.0000000 ## chi square test p-value for the null vs model
with(Bleaching_2, (null.deviance - deviance))## [1] 226.4038 with(Bleaching_2, (df.null - df.residual))## [1] 16 with(Bleaching_2, pchisq(null.deviance - deviance, df.null - df.residual, lower.tail = FALSE))## [1] 3.459688e-39 logLik(Bleaching_2)## 'log Lik.' -402.7962 (df=17) # model as a whole fits significantly better than an empty model.
aod::wald.test(b = coef(Bleaching_2), Sigma = vcov(Bleaching_2), Terms = 5:6)## Wald test:
## ----------
##
## Chi-squared test:
## X2 = 10.9, df = 2, P(> X2) = 0.0043 ## odds ratios and 95% CI
exp(cbind(OR = coef(Bleaching_2), confint(Bleaching_2)))## OR 2.5 % 97.5 %
## (Intercept) 0.8213709 0.27353838 2.2956273
## Z..m. 0.8789055 0.80442673 0.9594841
## Spp2GAPL 3.2728684 0.68787138 17.6540535
## Spp2PACL 2.3879275 0.68671983 8.8117814
## Spp2PAGI 0.5090179 0.13439142 2.0128705
## Spp2PAVA 3.0043606 0.87203319 11.0285384
## Spp2POLO 1.3027066 0.30814858 5.8181114
## Spp2POPA 1.6533905 0.44584873 6.4254639
## Spp2POCSP 0.2328654 0.06877018 0.8209344
## Year_F2016 0.2909686 0.19799682 0.4201602
## Z..m.:Spp2GAPL 0.9008468 0.69385639 1.1082381
## Z..m.:Spp2PACL 0.9820201 0.85790407 1.1127702
## Z..m.:Spp2PAGI 1.0084365 0.83827627 1.1593465
## Z..m.:Spp2PAVA 0.9222729 0.79363964 1.0488723
## Z..m.:Spp2POLO 0.8286930 0.62785800 1.0386842
## Z..m.:Spp2POPA 1.0335397 0.92602240 1.1510015
## Z..m.:Spp2POCSP 1.1023881 0.98433203 1.2335566 # Fitted values and plot Change in DProportion ~ days for different Intial D_Proportions
# create a new data frame
NewData <- expand.grid(Z..m.=0:28,
Year_F=unique(Bleaching_WP_SPP_2015_16$Year_F),
Spp2=unique(Bleaching_WP_SPP_2015_16$Spp2))
# get prediction
pred <- predict(Bleaching_2,newdata= NewData,se.fit=TRUE,
type = "response")
# CI function
make_ci <- function(pred, data){
# fit, lower, and upper CI
fit <- pred$fit
se<-pred$se.fit
lower <- fit - 1.96*pred$se.fit
upper <- fit + 1.96*pred$se.fit
return(data.frame(fit, lower, upper, data))
}
my_pred <- make_ci(pred, NewData)
Bleach_prop <- ggplot(my_pred, aes (Z..m., fit, colour=Spp2)) +
geom_line() + theme_bw( ) +
geom_ribbon(aes(ymin=my_pred$lower,
ymax=my_pred$upper),
linetype=0, alpha=0.1) +
scale_y_continuous(name=expression(Proportion~bleached),
expand = c(0, 0)) +
scale_x_continuous(name=expression("Depth (m)"),
breaks = seq(0,84,8),
expand = c(0, 0))
Bleach_prop+ facet_grid(~Year_F)
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.
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.
Grolemund, Garrett, and Hadley Wickham. 2011. “Dates and Times Made Easy with lubridate.” Journal of Statistical Software 40 (3): 1–25. https://www.jstatsoft.org/v40/i03/.
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/.
Spinu, Vitalie, Garrett Grolemund, and Hadley Wickham. 2021. Lubridate: Make Dealing with Dates a Little Easier. https://CRAN.R-project.org/package=lubridate.
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.
———. 2021. 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, 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.