###3_Martinez2019_coral_extension ###Martinez et al. 2019 ###www.cookbook-r.com library(ggplot2) library(plyr) coral <- read.csv("MartinezEtAl2019_CaribbeanCoralCalcification.csv") #select corals that survived experiment skeleton <- subset(coral, coral$status.survival == 0) #calculate extension rate (cm/year). Convert mm to cm. Convert 2-year extension (24 months) to 1 year skeleton$extension.rate_cm.year <- ((skeleton$extension_mm/10)/2) #make subsets by coral species p.ast <-subset (skeleton, skeleton$species == "Porites astreoides") sid <- subset(skeleton, skeleton$species == "Siderastrea") p.por <-subset(skeleton, skeleton$species =="Porites porites") pd <- position_dodge(0.5) #dodge overlapping objects 0.5 to the left and right #summarySE from http://www.cookbook-r.com/Manipulating_data/Summarizing_data/ #MODIFIED to include median and to exclude ci ## Summarizes data. ## Gives count, mean, median, standard deviation and standard error of the mean. ## data: a data frame. ## measurevar: the name of a column that contains the variable to be summariezed ## groupvars: a vector containing names of columns that contain grouping variables ## na.rm: a boolean that indicates whether to ignore NA's summarySE <- function(data=NULL, measurevar, groupvars=NULL, na.rm=FALSE, .drop=TRUE) { # New version of length which can handle NA's: if na.rm==T, don't count them length2 <- function (x, na.rm=FALSE) { if (na.rm) sum(!is.na(x)) else length(x) } # This does the summary. For each group's data frame, return a vector with # N, mean, median and sd datac <- ddply(data, groupvars, .drop=.drop, .fun = function(xx, col) { c(N = length2(xx[[col]], na.rm=na.rm), mean = mean (xx[[col]], na.rm=na.rm), median = median (xx[[col]], na.rm=na.rm), sd = sd (xx[[col]], na.rm=na.rm) ) }, measurevar ) # Rename the "mean" column datac <- rename(datac, c("mean" = measurevar)) datac$se <- datac$sd / sqrt(datac$N) # Calculate standard error of the mean return(datac) } ###P. astreoides hist(p.ast$extension.rate_cm.year) by(p.ast$extension.rate_cm.year, p.ast$origin, hist) #plot histogram for each origin site by(p.ast$extension.rate_cm.year, p.ast$transplant, hist) #plot histogram for each transplant site qqnorm(p.ast$extension.rate_cm.year) qqline(p.ast$extension.rate_cm.year) shapiro.test(p.ast$extension.rate_cm.year) #stats by origin. wilcox.test(extension.rate_cm.year ~ origin, data = p.ast) ext.p.ast.o.s <- summarySE(p.ast, measurevar ="extension.rate_cm.year", groupvars = c("species", "origin"), na.rm =T) ext.p.ast.o.s #stats by transplantation. wilcox.test(extension.rate_cm.year ~ transplant, data = p.ast) ext.p.ast.t.s <- summarySE(p.ast, measurevar ="extension.rate_cm.year", groupvars = c("species", "transplant"), na.rm =T) ext.p.ast.t.s #plot. ext.p.ast.or.t.s <- summarySE(p.ast, measurevar ="extension.rate_cm.year", groupvars = c("species", "origin", "transplant"), na.rm =T) ext.p.ast.or.t.s ext.p.ast.or.t.p <- ggplot(ext.p.ast.or.t.s, aes(x=transplant, y=extension.rate_cm.year, shape= origin)) + geom_errorbar(aes(ymin=extension.rate_cm.year-se, ymax=extension.rate_cm.year+se), width=.1, position = pd) + geom_line(position = pd) + geom_point(position = pd, size =3)+ xlab("Transplantation site") + ylab(expression(paste("Linear extension (cm ", year^{-1}, ")")))+ scale_shape_manual(values = c(16,1))+ ggtitle(expression(bolditalic("Porites astreoides")))+ theme(plot.title = element_text(size = 15), axis.text=element_text(size=12), axis.title=element_text(size=12), panel.background = element_blank(), panel.border = element_rect(fill = NA), legend.background = element_rect(linetype = "solid", colour = "black"), legend.key = element_rect(fill = "white")) ext.p.ast.or.t.p ###S.siderea hist(sid$extension.rate_cm.year) by(sid$extension.rate_cm.year, sid$origin, hist) #plot histogram for each origin site by(sid$extension.rate_cm.year, sid$transplant, hist) #plot histogram for each transplant site qqnorm(sid$extension.rate_cm.year) qqline(sid$extension.rate_cm.year) shapiro.test(sid$extension.rate_cm.year) #stats by origin. ext.sid.or.s <-summarySE(sid, measurevar = "extension.rate_cm.year", groupvars = c("species","origin"), na.rm = T) ext.sid.or.s kruskal.test(sid$extension.rate_cm.year ~ origin, data = sid) pairwise.wilcox.test(sid$extension.rate_cm.year, sid$origin) #stats by transplantation. ext.sid.t.s <- summarySE(sid, measurevar = "extension.rate_cm.year", groupvars = c("species","transplant"), na.rm = T) ext.sid.t.s wilcox.test(extension.rate_cm.year ~ transplant, data = sid) #plot. ext.sid.or.t.s <- summarySE(sid, measurevar = "extension.rate_cm.year", groupvars = c("species","origin", "transplant"), na.rm = T) ext.sid.or.t.s ext.sid.or.t.p <- ggplot(ext.sid.or.t.s, aes(x=transplant, y=extension.rate_cm.year, shape= origin)) + geom_errorbar(aes(ymin=extension.rate_cm.year-se, ymax=extension.rate_cm.year+se), width=.1, position = pd) + geom_line(position = pd) + geom_point(position = pd, size =3)+ xlab("Transplantation site") + ylab(expression(paste("Linear extension (cm ", year^{-1}, ")")))+ scale_shape_manual(values = c(16,1,17,16,1,17))+ ggtitle(expression(bolditalic("Siderastrea siderea"))) + theme(plot.title = element_text(size = 15), axis.text=element_text(size=12), axis.title=element_text(size=12), panel.background = element_blank(), panel.border = element_rect(fill = NA), legend.background = element_rect(linetype = "solid", colour = "black"), legend.key = element_rect(fill = "white")) ext.sid.or.t.p ###P.porites hist(p.por$extension.rate_cm.year) by(p.por$extension.rate_cm.year, p.por$transplant, hist) #plot histogram for each transplant site qqnorm(p.por$extension.rate_cm.year) qqline(p.por$extension.rate_cm.year) shapiro.test(p.por$extension.rate_cm.year) #stats by transplantation. wilcox.test(extension.rate_cm.year ~ transplant, data = p.por) #plot. ext.p.por.or.t.s <- summarySE(p.por, measurevar ="extension.rate_cm.year", groupvars = c("species", "origin", "transplant"), na.rm =T) ext.p.por.or.t.s ext.p.por.or.t.p <-ggplot(ext.p.por.or.t.s, aes(x=transplant, y=extension.rate_cm.year, shape= origin)) + geom_errorbar(aes(ymin=extension.rate_cm.year-se, ymax=extension.rate_cm.year+se), width=.1, position = pd) + geom_line(position = pd) + geom_point(position = pd, size =3)+ xlab("Transplantation site") + ylab(expression(paste("Linear extension (cm ", year^{-1}, ")")))+ scale_shape_manual(values = c(17))+ scale_y_continuous(labels = scales::number_format(accuracy = 0.1))+ ggtitle(expression(bolditalic("Porites porites"))) + theme(plot.title = element_text(size = 15), axis.text=element_text(size=12), axis.title=element_text(size=12), panel.background = element_blank(), panel.border = element_rect(fill = NA), legend.background = element_rect(linetype = "solid", colour = "black"), legend.key = element_rect(fill = "white")) ext.p.por.or.t.p