## R code to produce Figure 3, and the accompanying analyses presented in the text, in: ## Burgess SC, Johnston EC, Wyatt ASJ, Leichter JJ, Edmunds PJ (in press) Response diversity in corals: hidden differences in bleaching mortality among cryptic Pocillopora species. Ecology. ## Code written by Scott Burgess. Feb 2021. Send comments or corrections to sburgess@bio.fsu.edu ## R version 3.6.3 (2020-02-29) # Load required library library(plyr) # Import data # setwd("") dat <- read.csv("Data on Bleaching 2019.csv") # Prepare data dat$Site <- as.factor(dat$Site) dat$Quadrat <- as.factor(dat$Quadrat) dat$Site.Quadrat <- interaction(dat$Site,dat$Quadrat) dat$BleachAny <- ifelse(dat$Bleaching=="0",0,1) dat$Bleach1 <- ifelse(dat$Bleaching=="1",1, ifelse(dat$Bleaching=="0",0,NA)) dat$Bleach2 <- ifelse(dat$Bleaching=="2",1, ifelse(dat$Bleaching=="0",0,NA)) dat$Bleach3 <- ifelse(dat$Bleaching=="3",1, ifelse(dat$Bleaching=="0",0,NA)) # Subset data into each month datmarch <- dat[dat$Month=="March",] datmay <- dat[dat$Month=="May",] # Sample sizes # Size and bleaching were measured for 641 out of 1023 (63%) colonies assigned to a bleaching category in May 2019. length(datmay[!is.na(datmay$Longest.cm) & !is.na(datmay$BleachAny),which(names(datmay)=="Longest.cm")]) # 641 length(datmay$Bleaching[!is.na(datmay$Bleaching)]) # 1023 (641/1032) * 100 # 62.11% # Sample sizes at each site (from Figure 3 caption) # n=328 at Site 1 length(datmay[datmay$Site=="1" & !is.na(datmay$Longest.cm) & !is.na(datmay$BleachAny),which(names(datmay)=="Longest.cm")]) # n=313 at Site 2 length(datmay[datmay$Site=="2" & !is.na(datmay$Longest.cm) & !is.na(datmay$BleachAny),which(names(datmay)=="Longest.cm")]) # Prepare data for Figure 3 foo <- datmay[datmay$Site=="1",which(names(datmay)=="Longest.cm")] LengthvecS1 <- seq(min(foo,na.rm=T),max(foo,na.rm=T),length=100) foo <- datmay[datmay$Site=="2",which(names(datmay)=="Longest.cm")] LengthvecS2 <- seq(min(foo,na.rm=T),max(foo,na.rm=T),length=100) l <- c(LengthvecS1, LengthvecS2) newdat <- data.frame(Longest.cm=l,Site=as.factor(sort(rep(1:2,100)))) ## BleachAny in May m1 <- glm(BleachAny ~ Longest.cm * Site, data=datmay,family="binomial") p <- predict(m1,newdat,se.fit=T) pred.BleachAny.Length.May <- data.frame(newdat, mean=p$fit, upr = p$fit + 2*p$se.fit, lwr = p$fit - 2*p$se.fit) ## Bleach1 in May m1 <- glm(Bleach1 ~ Longest.cm * Site, data=datmay,family="binomial") p <- predict(m1,newdat,se.fit=T) pred.Bleach1.Length.May <- data.frame(newdat, mean=p$fit, upr = p$fit + 2*p$se.fit, lwr = p$fit - 2*p$se.fit) ## Bleach2 in May m1 <- glm(Bleach2 ~ Longest.cm * Site, data=datmay,family="binomial") p <- predict(m1,newdat,se.fit=T) pred.Bleach2.Length.May <- data.frame(newdat, mean=p$fit, upr = p$fit + 2*p$se.fit, lwr = p$fit - 2*p$se.fit) ## Bleach3 in May m1 <- glm(Bleach3 ~ Longest.cm * Site, data=datmay,family="binomial") p <- predict(m1,newdat,se.fit=T) pred.Bleach3.Length.May <- data.frame(newdat, mean=p$fit, upr = p$fit + 2*p$se.fit, lwr = p$fit - 2*p$se.fit) # Make Figure 3 labs <- c("Partial","Moderate","Severe","Any") x.ticks <- seq(0,45,5) panel.labs <- c("a) Site 1 (10m)","b) Site 2 (10m)") points.cex <- 1.5 ax.cex <- 1.2 x.ax <- seq(0,1,0.2) cols <- c("grey80","grey60","grey30","black") tcksize <- 0.05 xmin=0 xmax <- 45 brks <- seq(0,xmax,5) dev.new(width=8,height=4) nf <- layout(matrix(1:4,2,2,byrow=T),c(3,3),c(2,0.6,2,0.6)) # layout.show(n=4) par(oma=c(2,2,1,0),mar=c(0,0,0,0)) ## Site 1 par(mar=c(1,3,1,1)) plot(c(min(x.ticks),max(x.ticks)),c(0,1),type="n",ylab="",xlab="",yaxt="n",xaxt="n",bty="l") axis(side=1,at=x.ticks,cex.axis=ax.cex) axis(side=2,at=x.ax,las=1,cex.axis=ax.cex) mtext(panel.labs[1],side=3,line=-1,cex=ax.cex,adj=0.05) # Site 1 Bleach1 maxsize <- max(datmay[datmay$Site=="1" & !is.na(datmay$Longest.cm) & !is.na(datmay$Bleach1),which(names(datmay)=="Longest.cm")]) foo <- pred.Bleach1.Length.May[pred.Bleach1.Length.May$Site=="1" & pred.Bleach1.Length.May$Longest.cm