# Jeff Shima # Latest update: 25 November 2020 # This script generates analyses and figures that appear in 'Lunar rhythms in growth of larval fish '(Shima et al. doi:10.1098/rspb.2020.2609) # The script requires data available from Dryad: doi:10.5061/dryad.zgmsbcc93 # User intervention is also required at lines 36, 610, and 840 to specify working directory information ########################### # package libraries loaded library(dplyr) library(ggplot2) library(ggfortify) library(gridExtra) library(readr) library(lubridate) library(quantreg) library(tidyr) library(ggthemes) library(scales) library(car) library(igraph) library(ggExtra) library(cowplot) library(lme4) library(MuMIn) library(performance) library(effects) library(nlme) library(RColorBrewer) library(ggforce) rm(list=ls()) ########################### # file.choose() #opens a dialog box to locate/copy path to data (that can be pasted into quotes below) setwd("D:\\Dropbox\\Working Hard Drive MIRROR\\Research\\Papers in progress\\1 - Shima et al (PROC B)\\dryad data upload\\testing\\doi_10.5061_dryad.zgmsbcc93__v2") ########################### ############################### #SUMMARY OF APPROACH # (1) Constrain full data set to new settlers, to avoid confounding effects of post-settlement selection (and challenges associated with back-calculating size using Modified-Fry) # (2) Back-calculate fish size (TL, in mm) using the Modified Fry method, as per Wilson, Vigliola and Meeken (2009) # (3) Fit candidate growth models, choose best using AIC, obtain residuals as an estimate of de-trended (residual) size # (4) Calculate "residual growth" as change in residual size at each time step, for each fish # (5) Test for lunar periodicity in residual growth using periodic regression # (6) Explore potential explanations for lunar periodicity in residual growth using a secondary analysis ############################### # (1) Constrain full.data set to new settlers, to avoid confounding effects of post-settlement selection (and challenges associated with back-calculating size using Modified-Fry) # (2) Back-calculate fish size (TL, in mm) using the Modified Fry method, as per Wilson, Vigliola and Meeken (2009) ############################### # generating the datasets - step 1 ############################### full.data<-read.csv("thalassoma-hardwicke-otolith-growth-increments.csv") glimpse(full.data) # look at the data work.set <- filter(full.data,days.prior.to.set==0) %>% #constructing a data set containing Lcpt (length at capture) and Rcpt (otolith radius at capture); filter on days.prior.to.set gives an appropriate single row per fish; filter on stage constrains full data to new settlers filter(stage=="settler") %>% select (fish.id, stage, pld, birth.date, etoh.tl.est, total.otolith.radius, radius.width, date) %>% mutate(Lop = 1.5, #setting estimated larval length at hatch as 1.5 (based on literature estimates from related species) Rop = 2.5) %>% #setting a common otolith radius at hatch as 2.5 microns rename(Rcpt = total.otolith.radius) %>% rename(Lcpt = etoh.tl.est) %>% rename(Ri = radius.width) # otolith radius at settlement ############################### ############################### # implementing nonlinear regression to estimate parameters needed for MF back-calculation ############################### b.start <- 1.66 #starting value for estimation of b (these are values from a SAS nLin procedure) c.start <- 0.4 #starting value for estimation of c (these are values from a SAS nLin procedure) # model m.set <- nls(Lcpt ~ Lop-(b*(Rop**c))+(b*(Rcpt**c)), start=list(b = b.start, c = c.start), data=work.set) #estimating parameters using both settlers only #estimated parameters summary(m.set) #estimates are: b=1.74880, c=0.39147 ############################### ############################### # generating the datasets - step 2 ############################### work.set <- work.set %>% select(fish.id, Rcpt) %>% left_join(full.data, by= "fish.id") %>% mutate(Lop = 1.5, #setting estimated larval length at hatch as 1.5 (based on literature estimates from related species) Rop = 2.5) %>% #setting a common otolith radius at hatch as 2.5 microns rename(Lcpt = etoh.tl.est) %>% rename(Ri = radius.width) %>% # otolith radius at settlement mutate(b=1.74880, # from nls() fit above c=0.39147 , # from nls() fit above a=Lop-b*Rop^c, Li=a+exp( (log(Lop-a)) + ( (log(Lcpt-a)) - (log(Lop-a)) ) * ( (log(Ri)-log(Rop) ) / ( log(Rcpt)-log(Rop) ) ) ) , Pred.L.cpt= a +b*(Rcpt^c), lag.Li = lag(Li, 1), delta.Li = if_else(Li-lag.Li>0,Li-lag.Li,NA_real_)) %>% #calculating the daily growth rate (in mm TL), using the Modified-Fry Method (Wilson and Vigliola paper) select (fish.id, increment, incr.width, Ri, larval.age, pld, Li, lag.Li, delta.Li, date, lunar.day, lunar.month,lunar.phase.category.centered, illumination, cloud.cover.airs.nighttime.3x3, time.moonrise, time.moonset, time.moonrise.1, time.meridian.passing) # include additional environmental data as required ############################### env.spawning <- filter(full.data,increment==1) %>% #constructing descriptors of birth date filter(stage=="settler") %>% select(fish.id, lunar.day, lunar.month, lunar.day.bin, lunar.phase.category.start, lunar.phase.category.centered) %>% rename(LD.birth=lunar.day, LD.birth.bin=lunar.day.bin, LM.birth=lunar.month, lunar.phase.birth.start=lunar.phase.category.start, lunar.phase.birth.centered=lunar.phase.category.centered,) %>% mutate(LD.birth.squared = LD.birth**2) ############################### work.set <- left_join(work.set, env.spawning, by="fish.id") %>% filter(!is.na(lunar.phase.birth.start)) %>% #removing observations for which we cannot align to a lunar calendar (e.g., unable to age post-settlement increments) mutate(Li.2=if_else(larval.age==pld,NA_real_,Li)) %>% # created a set of Li that excludes length on final day of larval development. This will used to generate residual size that at age x for only the set of individuals present at x+1 (to facilitate estimates of delta.resid that are not biased by settlement processes) mutate(date=ymd(date)) %>% #correcting date format mutate(birth.date=ymd(date)) %>% #correcting date format mutate(set.date=ymd(date)) %>% #correcting date format mutate(time.moonrise = hm(time.moonrise)) %>% #correcting time format mutate(time.moonset = hm(time.moonset)) %>% #correcting time format mutate(time.moonrise.1 = hm(time.moonrise.1)) %>% #correcting time format mutate(time.meridian.passing = hms(time.meridian.passing)) glimpse(work.set) # look at the data ############################### ############################### ############################### # (3) Fit candidate growth models, choose best using AIC, obtain residuals as an estimate of de-trended (residual) size ############################### ############################### # VBGF k.start <- 0.05 #starting value for estimation Linf.start <- 15 #starting value for estimation a0.start <- 1 #starting value for estimation m.vbgf <- nls(Li ~ Linf*(1-exp(-k*(larval.age-a0))), start=list(Linf = Linf.start, k = k.start, a0 = a0.start), data=work.set) summary(m.vbgf) # Gompertz a.start <- 15 b.start <- 0.05 c.start <- -0.2 m.gompertz <- nls(Li ~ a*(exp(-b*(exp(-c*larval.age)))), start=list(a = a.start, b = b.start, c = c.start), data=work.set) summary(m.gompertz) # Logistic L.start <- 15 x0.start <- 15 k.start <- .5 m.logistic <- nls(Li ~ L/(1+exp(-k*(larval.age-x0))), start = list(L=L.start, x0=x0.start, k=k.start), data=work.set) summary(m.logistic) # Quadratic b1.start <- -.5 b2.start <- .3 b3.start <- -1 m.quadratic <- nls(Li ~ b1*larval.age**2 + b2*larval.age + b3, start = list(b1=b1.start, b2=b2.start, b3=b3.start), data=work.set) summary(m.quadratic) # Unstructured m.larval.age <- lm(Li ~ factor(larval.age), data = work.set) #fitting model with larval age as a factor (this effectively centers Li for each larval age) summary(m.larval.age) ############################### #Identifying best model AICc(m.vbgf, m.gompertz, m.logistic, m.quadratic, m.larval.age) #Unconstrained model has lowest AIC # fitting the best model to Li.2 (which excludes size on final day of larval life), to facilitate unbiased estimates of resid.growth m.larval.age.2 <- lm(Li.2 ~ factor(larval.age), data = work.set, na.action = na.exclude) #using "na.action = na.exclude" argument here to pad the object with NAs in the correct positions to have the same number of rows as the original data frame. summary(m.larval.age.2) #Appending residuals from top two models vbgf.resid <- resid(m.vbgf) work.set$resid <- vbgf.resid larval.age.resid <- resid(m.larval.age) work.set$larval.age.resid <- larval.age.resid glimpse(work.set) larval.age.predict <- predict(m.larval.age) work.set$larval.age.predict <- larval.age.predict glimpse(work.set) larval.age.resid.2 <- resid(m.larval.age.2) work.set$larval.age.resid.2 <- larval.age.resid.2 glimpse(work.set) larval.age.predict.2 <- predict(m.larval.age.2) work.set$larval.age.predict.2 <- larval.age.predict.2 glimpse(work.set) ############################### ############################### # Plotting the best model for de-trending size-at-age (and other functions, for appendix figure) ############################### ############################### #Fig A1 ############################### min.age<-min(work.set$larval.age) #using these functions to delineate the range of the new x variable for fitted line max.age<-max(work.set$larval.age) new.x<- expand.grid(larval.age = seq(min.age,max.age, length=1000)) new.y = predict(m.vbgf, newdata = new.x, se.fit = TRUE) new.y = data.frame(new.y) head(new.y) addThese<- data.frame(new.x,new.y) %>% #values to facilitate fitted VBGF on plot rename(Li = new.y) fig.larval.age.growth <- ggplot(work.set, aes(x = larval.age, y = Li)) + geom_point(work.set, mapping = aes(x = larval.age, y = Li), size = .3, alpha = 0.3, show.legend=FALSE, position = position_jitter(width=.3) ) + geom_smooth(data=addThese, stat='identity', size=1.1,colour="black") + #adding VBGF geom_smooth(aes(x = larval.age, y = larval.age.predict), stat='identity', size=1.1,colour="red") + ylab("Larval size (TL, mm)") + xlab("Larval age (days)") + #relabeling x and y axes scale_x_continuous(limits = c(0,60), breaks = seq(from = 0, to = 60, by = 10)) + scale_y_continuous(limits = c(0,16), breaks = seq(from=0, to=15, by=3)) + theme_classic (base_size = 28) + theme(plot.title = element_text(colour="black", size=rel(1), hjust=0.05, vjust=-12)) fig.larval.age.growth # ggsave(plot = fig.larval.age.growth, width=10, height=8, dpi=300,filename = "figures\\Fig A1.pdf") # ggsave(plot = fig.larval.age.growth, width=10, height=8, dpi=300,filename = "figures\\Fig A1.png") #predicted size-at-settlement size.at.settlement <- work.set %>% #generating estimates to present in results section (to compare with previously reported empirical estimates of size-at-settlement) group_by(fish.id)%>% summarise(set.size=max(Li))%>% summarise(mean.set.size=mean(set.size), min.set.size=min(set.size), max.set.size=max(set.size), sd.set.size=sd(set.size)) ############################### #Appendix Fig A2 ############################### min.age<-min(work.set$larval.age) #using these functions to delineate the range of the new x variable for fitted line max.age<-max(work.set$larval.age) new.x<- expand.grid(larval.age = seq(min.age,max.age, length=1000)) new.y = predict(m.vbgf, newdata = new.x, se.fit = TRUE) new.y = data.frame(new.y) head(new.y) addThese<- data.frame(new.x,new.y) %>% rename(Li = new.y) figA2.vbgf.growth <- ggplot(work.set, aes(x = larval.age, y = Li)) + geom_point(work.set, mapping = aes(x = larval.age, y = Li), size = .3, alpha = 0.3, show.legend=FALSE, position = position_jitter(width=.3) ) + geom_smooth(data=addThese, stat='identity', size=1.1,colour="red") + ylab("Larval size (TL, mm)") + xlab("Larval age (days)") + ggtitle("(a) von Bertalanffy growth function (AIC: 24821)")+ scale_x_continuous(limits = c(0,60), breaks = seq(from = 0, to = 60, by = 10)) + scale_y_continuous(limits = c(0,16), breaks = seq(from=0, to=15, by=3)) + theme_classic (base_size = 16) + theme(plot.title = element_text(colour="black", size=rel(1), hjust=0.13, vjust=-9)) figA2.vbgf.growth new.x<- expand.grid(larval.age = seq(min.age,max.age, length=1000)) new.y = predict(m.logistic, newdata = new.x, se.fit = TRUE) new.y = data.frame(new.y) head(new.y) addThese<- data.frame(new.x,new.y) %>% rename(Li = new.y) figA2.logistic.growth <- ggplot(work.set, aes(x = larval.age, y = Li)) + geom_point(work.set, mapping = aes(x = larval.age, y = Li), size = .3, alpha = 0.3, show.legend=FALSE, position = position_jitter(width=.3) ) + geom_smooth(data=addThese, stat='identity', size=1.1,colour="red") + ylab("Larval size (TL, mm)") + xlab("Larval age (days)") + ggtitle("(b) Logistic growth function (AIC: 31040)")+ scale_x_continuous(limits = c(0,60), breaks = seq(from = 0, to = 60, by = 10)) + scale_y_continuous(limits = c(0,16), breaks = seq(from=0, to=15, by=3)) + theme_classic (base_size = 16) + theme(plot.title = element_text(colour="black", size=rel(1), hjust=0.13, vjust=-9)) figA2.logistic.growth new.x<- expand.grid(larval.age = seq(min.age,max.age, length=1000)) new.y = predict(m.gompertz, newdata = new.x, se.fit = TRUE) new.y = data.frame(new.y) head(new.y) addThese<- data.frame(new.x,new.y) %>% rename(Li = new.y) figA2.gompertz.growth <- ggplot(work.set, aes(x = larval.age, y = Li)) + geom_point(work.set, mapping = aes(x = larval.age, y = Li), size = .3, alpha = 0.3, show.legend=FALSE, position = position_jitter(width=.3) ) + geom_smooth(data=addThese, stat='identity', size=1.1,colour="red") + ylab("Larval size (TL, mm)") + xlab("Larval age (days)") + ggtitle("(c) Gompertz growth function (AIC: 27922)")+ scale_x_continuous(limits = c(0,60), breaks = seq(from = 0, to = 60, by = 10)) + scale_y_continuous(limits = c(0,16), breaks = seq(from=0, to=15, by=3)) + theme_classic (base_size = 16) + theme(plot.title = element_text(colour="black", size=rel(1), hjust=0.13, vjust=-9)) figA2.gompertz.growth new.x<- expand.grid(larval.age = seq(min.age,max.age, length=1000)) new.y = predict(m.quadratic, newdata = new.x, se.fit = TRUE) new.y = data.frame(new.y) head(new.y) addThese<- data.frame(new.x,new.y) %>% rename(Li = new.y) figA2.quadratic.growth <- ggplot(work.set, aes(x = larval.age, y = Li)) + geom_point(work.set, mapping = aes(x = larval.age, y = Li), size = .3, alpha = 0.3, show.legend=FALSE, position = position_jitter(width=.3) ) + geom_smooth(data=addThese, stat='identity', size=1.1,colour="red") + ylab("Larval size (TL, mm)") + xlab("Larval age (days)") + ggtitle("(d) Quadratic growth function (AIC: 26875)")+ scale_x_continuous(limits = c(0,60), breaks = seq(from = 0, to = 60, by = 10)) + scale_y_continuous(limits = c(0,16), breaks = seq(from=0, to=15, by=3)) + theme_classic (base_size = 16) + theme(plot.title = element_text(colour="black", size=rel(1), hjust=0.13, vjust=-9)) figA2.quadratic.growth figA2.larval.age.growth <- ggplot(work.set, aes(x = larval.age, y = Li)) + geom_point(work.set, mapping = aes(x = larval.age, y = Li), size = .3, alpha = 0.3, show.legend=FALSE, position = position_jitter(width=.3) ) + geom_smooth(aes(x = larval.age, y = larval.age.predict), stat='identity', size=1.1,colour="red") + ylab("Larval size (TL, mm)") + xlab("Larval age (days)") + ggtitle(" (e) Unrestricted age-dependent growth (AIC: 23537)")+ scale_x_continuous(limits = c(0,60), breaks = seq(from = 0, to = 60, by = 10)) + scale_y_continuous(limits = c(0,16), breaks = seq(from=0, to=15, by=3)) + theme_classic (base_size = 16) + theme(plot.title = element_text(colour="black", size=rel(1), hjust=0.13, vjust=-9)) figA2.larval.age.growth figA2.multipanel<- grid.arrange(figA2.vbgf.growth, figA2.logistic.growth, figA2.gompertz.growth, figA2.quadratic.growth, figA2.larval.age.growth, ncol = 2, nrow = 3) # ggsave(plot = figA2.multipanel, width=12, height=13.5, dpi=300,filename = "figures\\Fig A2.pdf") # ggsave(plot = figA2.multipanel, width=12, height=13.5, dpi=300,filename = "figures\\Fig A2.png") ############################### #Fig A3 ############################### fig.vbgf.residuals <- ggplot(work.set, aes(x = larval.age, y = resid)) + geom_point(work.set, mapping = aes(x = larval.age, y = resid), size = .3, alpha = 0.3, show.legend=FALSE, position = position_jitter(width=.3) ) + geom_smooth(stat="smooth",position="identity", method = "auto", level= 0.95, size=1.1,colour="black") + #for description of this fitted line, see: https://ggplot2.tidyverse.org/reference/geom_smooth.html geom_smooth(aes(y=0), stat='identity', size=.5, colour="black", linetype="dashed") + ylab("Residual larval size (TL, mm)") + xlab("Larval age (days)") + ggtitle(" (a) von Bertalanffy growth model")+ scale_x_continuous(limits = c(0,60), breaks = seq(from = 0, to = 60, by = 10)) + theme_classic (base_size = 28) + theme(plot.title = element_text(colour="black", size=rel(1), hjust=0.05, vjust=-5)) fig.vbgf.residuals fig.larval.age.resid <- ggplot(work.set, aes(x = larval.age, y = larval.age.resid)) + geom_point(work.set, mapping = aes(x = larval.age, y = larval.age.resid), size = .3, alpha = 0.3, show.legend=FALSE, position = position_jitter(width=.3) ) + geom_smooth(stat="smooth",position="identity", method = "auto", level= 0.95, size=1.1,colour="black") + #for description of this fitted line, see: https://ggplot2.tidyverse.org/reference/geom_smooth.html geom_smooth(aes(y=0), stat='identity', size=.5, colour="black", linetype="dashed") + ylab("Residual larval size (TL, mm)") + xlab("Larval age (days)") + ggtitle(" (b) Unrestricted age-dependent growth model")+ scale_x_continuous(limits = c(0,60), breaks = seq(from = 0, to = 60, by = 10)) + theme_classic (base_size = 28) + theme(plot.title = element_text(colour="black", size=rel(1), hjust=0.05, vjust=-5)) fig.larval.age.resid figA3.multipanel<- grid.arrange(fig.vbgf.residuals, fig.larval.age.resid, ncol = 2, nrow = 1) # ggsave(plot = figA3.multipanel, width=20, height=8, dpi=300,filename = "figures\\Fig A3.pdf") # ggsave(plot = figA3.multipanel, width=20, height=8, dpi=300,filename = "figures\\Fig A3.png") ############################### ############################### # (4) Calculate "residual growth" as change in residual size at each time step, for each fish ############################### work.set <- work.set %>% mutate(lag.larval.age.resid.2 = lag(larval.age.resid.2, 1), #calculating residual growth from "unrestricted" model (i.e., the best model) delta.resid = if_else(Li-lag.Li>0,larval.age.resid-lag.larval.age.resid.2,NA_real_)) #diagnostic plot of delta.resid to confirm that the use of "larval.age.resid.2" results in unbiased estimates of delta.resid (i.e., centered on 0) ggplot(work.set, aes(x = larval.age, y = delta.resid)) + geom_point(size = 2, alpha = .2, position = position_dodge(width=2)) + geom_smooth(work.set, mapping = aes(x = larval.age, y = delta.resid), stat="smooth",position="identity", span=10, method = "auto", level= 0.95, formula = y ~ x, size=1.1, colour="red") work.set1a <- work.set %>% filter(!is.na(delta.resid)) %>% group_by(lunar.phase.birth.centered, lunar.day) %>% summarise(mean.delta.resid = mean(delta.resid), n=n(), sd=sd(delta.resid), se=sd/sqrt(n)) work.set1b <- work.set %>% filter(!is.na(delta.resid)) %>% group_by(lunar.phase.birth.centered, larval.age) %>% summarise(mean.delta.resid = mean(delta.resid), n=n(), sd=sd(delta.resid), se=sd/sqrt(n)) work.set1c <- work.set %>% group_by(fish.id) %>% filter(larval.age==3) ############################################ #Fig 1 ############################################ fig1a <- ggplot(data = work.set1c, aes(x = LD.birth, fill= lunar.phase.birth.centered, alpha=1 )) + geom_histogram(binwidth=1, colour="black") + ylab("Frequency") + xlab("Lunar day of birth") + scale_colour_brewer(palette = "Dark2") + scale_fill_brewer(palette = "Dark2") + scale_y_continuous(limits = c(0,40), breaks = seq(from=0, to=40, by=10)) + coord_cartesian(ylim=c(0,40)) + ggtitle("(a) Lunar cohorts") + theme_classic (base_size = 28) + theme(plot.title = element_text(colour="black", size=rel(1), hjust=.05, vjust=-4), legend.position = "none", axis.title.x = element_text(margin = margin(t = 64, r = 0, b = 20, l = 0)), legend.title=element_blank(), plot.margin = margin(t = 0, r = 20, b = 0, l = 0)) fig1a fig1a <- ggdraw(fig1a) + draw_image("figures\\moonphase.png", x= 0.042, y = -.36, scale=.84) fig1a brewer.pal(n = 4, name = "Dark2") #requesting the colour codes for Brewer "Dark2" palette fig1b <- ggplot(work.set1b, aes(x = larval.age, y = mean.delta.resid, colour=lunar.phase.birth.centered)) + geom_point(size = 5, alpha = .4, position = position_dodge(width=2)) + geom_errorbar(aes(ymin = mean.delta.resid - se, ymax= mean.delta.resid+se),position = position_dodge(width=2), width = 0.5, size=1.1) + geom_smooth(work.set, mapping = aes(x = larval.age, y = delta.resid, colour=lunar.phase.birth.centered, fill=lunar.phase.birth.centered), stat="smooth",position="identity", span=3, method = "auto", level= 0.95, size=1.1, alpha=0.2) + geom_smooth(aes(y=0), stat='identity', size=1, colour="black", linetype="dashed") + ylab("Residual growth (mm/d)") + xlab("Larval age (days)\n") + ggtitle("(b) Growth trajectories of lunar cohorts")+ scale_colour_brewer(palette = "Dark2") + scale_fill_brewer(palette = "Dark2") + scale_x_continuous(limits = c(8,60), breaks = seq(from = 0, to = 60, by = 10)) + scale_y_continuous(limits = c(-0.1,0.1), breaks = seq(from=-.1, to=.1, by=0.05)) + annotate(geom = "point", x = 30, y = -0.06, colour = "#1B9E77", size=13, alpha=1) + annotate(geom = "point", x = 37, y = -0.06, colour = "#E7298A", size=13, alpha=1) + annotate(geom = "point", x = 44, y = -0.06, colour = "#7570B3", size=13, alpha=1) + annotate(geom = "point", x = 51, y = -0.06, colour = "#D95F02", size=13, alpha=1) + annotate(geom = "point", x = 58, y = -0.06, colour = "#1B9E77", size=13, alpha=1) + annotate(geom = "rect", xmin = 30, xmax =33, ymin =-0.05, ymax=-0.08, colour = "#FFFFFF", fill= "#FFFFFF", alpha=1) + annotate(geom = "rect", xmin = 37, xmax =40, ymin =-0.05, ymax=-0.08, colour = "#FFFFFF", fill= "#FFFFFF", alpha=1) + annotate(geom = "rect", xmin = 44, xmax =47, ymin =-0.05, ymax=-0.08, colour = "#FFFFFF", fill= "#FFFFFF", alpha=1) + annotate(geom = "rect", xmin = 51, xmax =52.5, ymin =-0.05, ymax=-0.08, colour = "#FFFFFF", fill= "#FFFFFF", alpha=1) + annotate(geom = "rect", xmin = 58, xmax =60, ymin =-0.05, ymax=-0.08, colour = "#FFFFFF", fill= "#FFFFFF", alpha=1) + annotate(geom = "curve", x = 29.9, xend=29.9, y=-.0520, yend=-.0680, curvature=-1, size=1, colour="#1B9E77") + annotate(geom = "curve", x = 36.9, xend=36.9, y=-.0520, yend=-.0680, curvature=-1, size=1, colour="#E7298A") + annotate(geom = "curve", x = 43.9, xend=43.9, y=-.0520, yend=-.0680, curvature=-1, size=1, colour="#7570B3") + annotate(geom = "curve", x = 50.9, xend=50.9, y=-.0520, yend=-.0680, curvature=-1, size=1, colour="#D95F02") + annotate(geom = "curve", x = 57.9, xend=57.9, y=-.0520, yend=-.0680, curvature=-1, size=1, colour="#1B9E77") + annotate(geom = "point", x = 23, y = 0.05, colour = "#E7298A", size=13, alpha=1) + annotate(geom = "point", x = 30, y = 0.05, colour = "#7570B3", size=13, alpha=1) + annotate(geom = "point", x = 37, y = 0.05, colour = "#D95F02", size=13, alpha=1) + annotate(geom = "point", x = 44, y = 0.05, colour = "#1B9E77", size=13, alpha=1) + annotate(geom = "point", x = 51, y = 0.05, colour = "#E7298A", size=13, alpha=1) + annotate(geom = "rect", xmin = 20, xmax =23, ymin =0.04, ymax=0.08, colour = "#FFFFFF", fill= "#FFFFFF", alpha=1) + annotate(geom = "rect", xmin = 27, xmax =30, ymin =0.04, ymax=0.08, colour = "#FFFFFF", fill= "#FFFFFF", alpha=1) + annotate(geom = "rect", xmin = 34, xmax =37, ymin =0.04, ymax=0.08, colour = "#FFFFFF", fill= "#FFFFFF", alpha=1) + annotate(geom = "rect", xmin = 41, xmax =44, ymin =0.04, ymax=0.08, colour = "#FFFFFF", fill= "#FFFFFF", alpha=1) + annotate(geom = "rect", xmin = 48, xmax =51, ymin =0.04, ymax=0.08, colour = "#FFFFFF", fill= "#FFFFFF", alpha=1) + annotate(geom = "curve", x = 23.1, xend=23.1, y=.0420, yend=.0580, curvature=-1, size=1, colour="#E7298A") + annotate(geom = "curve", x = 30.1, xend=30.1, y=.0420, yend=.0580, curvature=-1, size=1, colour="#7570B3") + annotate(geom = "curve", x = 37.1, xend=37.1, y=.0420, yend=.0580, curvature=-1, size=1, colour="#D95F02") + annotate(geom = "curve", x = 44.1, xend=44.1, y=.0420, yend=.0580, curvature=-1, size=1, colour="#1B9E77") + annotate(geom = "curve", x = 51.1, xend=51.1, y=.0420, yend=.0580, curvature=-1, size=1, colour="#E7298A") + annotate("text", x = 44, y = -0.082, label = "first quarter moons", size=7.5) + annotate("text", x = 37, y = 0.072, label = "last quarter moons", size=7.5) + theme_classic (base_size = 28) + theme(plot.title = element_text(colour="black", size=rel(1), hjust=.11, vjust=-7), legend.title=element_blank(), legend.background = element_rect(colour=NA, fill=NA), axis.title.x = element_text(margin = margin(t = 10, r = 0, b = 20, l = 0)), legend.position = c(0.17, 0.23), plot.margin = margin(t = 0, r = 0, b = 20, l = 10)) fig1b fig1 <- grid.arrange(fig1a, fig1b, ncol = 2, nrow = 1) #ggsave(plot = fig1, width=20, height=8, dpi=300,filename = "figures\\Fig 1.pdf") #ggsave(plot = fig1, width=20, height=8, dpi=300,filename = "figures\\Fig 1.png") ########################################################### ########################################################### # (5) Test for lunar periodicity in residual growth using periodic regression ########################################################### work.set2 <- work.set %>% filter(!is.na(delta.resid)) %>% filter(larval.age>=25) %>% #Constraining data to range of larval ages where (i)VBGF fits well, (2) variation in residual growth (Fig2b) is evident, and (3) sufficient samples exist to facilitate reasonable estimates (i.e., data are limited above age 50d) filter(larval.age<=50) %>% mutate(theta=(360*lunar.day)/28.5, #Adding parameters to facilitate periodic regression; 28.5d because numbering starts at 0 sin.theta=sin(theta*pi/180), cos.theta=cos(theta*pi/180)) # Approach to model selection as follows: # (1) use single (desired) fixed effects model structure, fit competing random effects structures using REML, and compare with AICc # (2) use winning random effect structure, fit competing fixed effects structures using maximum likelihood, compare with AICc # (3) refit best mixed-effect model using REML for unbiased parameter estimates m.lunar.rand1 <- lmer(delta.resid ~ sin.theta + (1|fish.id) , data = work.set2) m.lunar.rand2 <- lmer(delta.resid ~ sin.theta + (1+sin.theta|fish.id), data = work.set2) m.lunar.rand3 <- lmer(delta.resid ~ sin.theta + (1|date) , data = work.set2) m.lunar.rand4 <- lmer(delta.resid ~ sin.theta + (1+sin.theta|date) , data = work.set2) m.lunar.rand5 <- lmer(delta.resid ~ sin.theta + (1|fish.id) + (1|date) , data = work.set2) m.lunar.rand6 <- lmer(delta.resid ~ sin.theta + (1+sin.theta|fish.id) + (1|date) , data = work.set2) m.lunar.rand7 <- lmer(delta.resid ~ sin.theta + (1|fish.id) + (1+sin.theta|date) , data = work.set2) m.lunar.rand8 <- lmer(delta.resid ~ sin.theta + (1+sin.theta|fish.id) + (1+sin.theta|date) , data = work.set2) #the best random effects structure AICc(m.lunar.rand1,m.lunar.rand2,m.lunar.rand3,m.lunar.rand4,m.lunar.rand5,m.lunar.rand6,m.lunar.rand7,m.lunar.rand8) m.lunar.fixed1 <- lmer(delta.resid ~ sin.theta + (1+sin.theta|fish.id) + (1+sin.theta|date) , data = work.set2, REML=FALSE) m.lunar.fixed2 <- lmer(delta.resid ~ sin.theta + cos.theta + (1+sin.theta|fish.id) + (1+sin.theta|date) , data = work.set2, REML=FALSE) AICc(m.lunar.fixed1,m.lunar.fixed2) #the two models are essentially equivalent, will select the simpler model. m.lunar.final <- lmer(delta.resid ~ sin.theta + (1+sin.theta|fish.id) + (1+sin.theta|date) , data = work.set2) plot(m.lunar.final) qqnorm(resid(m.lunar.final)) qqline(resid(m.lunar.final)) anova(m.lunar.final) summary(m.lunar.final) #Values to report m.lunar <- lm(delta.resid ~ sin.theta , data = work.set2) anova(m.lunar) summary(m.lunar) ########################################################### #Figure 2 #Figure 2b min.day<-min(work.set2$lunar.day) max.day<-max(work.set2$lunar.day) new.x<- expand.grid(lunar.day = seq(min.day,max.day, length=1000)) new.x<- mutate(new.x,theta=(360*lunar.day)/28.5, sin.theta=sin(theta*pi/180), cos.theta=cos(theta*pi/180), mean.delta.resid=-0.0007310+-0.0129640*sin.theta) work.set2a <- work.set2 %>% group_by(lunar.day) %>% summarise(mean.delta.resid = mean(delta.resid), n=n(), sd=sd(delta.resid), se=sd/sqrt(n), conf=1.96*se) fig2b<-ggplot(work.set2a, aes(x = lunar.day, y = mean.delta.resid)) + geom_point(size = 5, alpha = .8) + geom_errorbar(aes(ymin = mean.delta.resid - se, ymax= mean.delta.resid+se), width = 0.5, size=1.1) + geom_smooth(data=new.x, aes(x = lunar.day, y = mean.delta.resid), stat='identity', size=1.1,colour="black") + ylab("Residual growth (mm/d)\n") + xlab("Lunar day") + ggtitle(" (b) Age-independent larval growth")+ scale_x_continuous(limits = c(-1,30), breaks = seq(from = 0, to = 30, by = 5)) + theme_classic (base_size = 27) + theme(axis.title.y = element_text(margin = margin(t = 0, r = 28, b = 0, l = 0)), plot.title = element_text(colour="black", size=rel(1), hjust=0.05, vjust=-3), axis.title.x = element_text(margin = margin(t = 55, r = 0, b = 0, l = 0))) fig2b <- ggdraw(fig2b) + draw_image("figures\\moonphase.png", x= 0.098, y = -.38, scale=.72) fig2b #Figure 2a moonlight<-read.csv("moonlight_ancillary_fig_2a.csv") glimpse(moonlight) #excel file of same name shows how this file was constructed fig2a <- ggplot(moonlight, mapping = aes(x = lunar.day, y = order)) + geom_raster(aes(fill=illum), stat="identity", hjust=0, vjust=0)+ scale_fill_gradient(name="", low="black", high="white")+ ylab("Portion of night illuminated") + xlab("Lunar day") + ggtitle(" (a) Nocturnal brightness across the lunar cycle")+ scale_x_continuous(limits = c(0,30), breaks = seq(from = 0, to = 30, by = 5)) + scale_y_continuous(limits = c(0,12), breaks = seq(from=0, to=12, by=3), labels=c("dawn", "03:00", "midnight", "21:00", "dusk")) + coord_cartesian(ylim=c(0,12)) + theme_classic (base_size = 27) + theme(plot.title = element_text(colour="black", size=rel(1), hjust=.05, vjust=-3), legend.position = 'none', axis.title.y = element_text(margin = margin(t = 0, r = 10, b = 0, l = 0))) fig2a fig2<- grid.arrange(fig2a, fig2b, ncol = 1, nrow = 2) # ggsave(plot = fig2, width=10, height=16, dpi=300,filename = "figures\\Fig 2.pdf") # ggsave(plot = fig2, width=10, height=16, dpi=300,filename = "figures\\Fig 2.png") ########################################################### ########################################################### #Appendix Fig A4 ########################################################### work.set3a <- work.set %>% filter(!is.na(delta.resid)) %>% filter(larval.age>=0) %>% filter(larval.age<=100) work.set3a1 <- work.set3a %>% group_by(lunar.day) %>% summarise(mean.delta.resid = mean(delta.resid), n=n(), sd=sd(delta.resid), se=sd/sqrt(n)) figA4a <- ggplot(work.set3a1, aes(x = lunar.day, y = mean.delta.resid)) + geom_point(size = 5, alpha = .8) + geom_errorbar(aes(ymin = mean.delta.resid - se, ymax= mean.delta.resid+se), width = 0.5, size=1.1) + geom_smooth(work.set3a, mapping = aes(x = lunar.day, y = delta.resid), stat="smooth",position="identity", method = "auto", level= 0.95, size=1.1, colour="black") + ylab("Residual growth (mm/d)\n") + xlab("Lunar day") + ggtitle(" (a) All larval ages")+ scale_x_continuous(limits = c(-1,30), breaks = seq(from = 0, to = 30, by = 5)) + theme_classic (base_size = 27) + theme(axis.title.y = element_text(margin = margin(t = 0, r = 28, b = 0, l = 0)), plot.title = element_text(colour="black", size=rel(1), hjust=0.05, vjust=-3)) figA4a work.set3b <- work.set %>% filter(!is.na(delta.resid)) %>% filter(larval.age>=25) %>% filter(larval.age<=50) work.set3b1 <- work.set3b %>% group_by(lunar.day) %>% summarise(mean.delta.resid = mean(delta.resid), n=n(), sd=sd(delta.resid), se=sd/sqrt(n)) figA4b <- ggplot(work.set3b1, aes(x = lunar.day, y = mean.delta.resid)) + geom_point(size = 5, alpha = .8) + geom_errorbar(aes(ymin = mean.delta.resid - se, ymax= mean.delta.resid+se), width = 0.5, size=1.1) + geom_smooth(work.set3b, mapping = aes(x = lunar.day, y = delta.resid), stat="smooth",position="identity", method = "auto", level= 0.95, size=1.1, colour="black") + ylab("Residual growth (mm/d)\n") + xlab("Lunar day") + ggtitle(" (b) Larval ages: 25-50d")+ scale_x_continuous(limits = c(-1,30), breaks = seq(from = 0, to = 30, by = 5)) + theme_classic (base_size = 27) + theme(axis.title.y = element_text(margin = margin(t = 0, r = 28, b = 0, l = 0)), plot.title = element_text(colour="black", size=rel(1), hjust=0.05, vjust=-3)) figA4b work.set3c <- work.set %>% filter(!is.na(delta.resid)) %>% filter(larval.age>=3) %>% filter(larval.age<=10) work.set3c1 <- work.set3c %>% group_by(lunar.day) %>% summarise(mean.delta.resid = mean(delta.resid), n=n(), sd=sd(delta.resid), se=sd/sqrt(n)) figA4c <- ggplot(work.set3c1, aes(x = lunar.day, y = mean.delta.resid)) + geom_point(size = 5, alpha = .8) + geom_errorbar(aes(ymin = mean.delta.resid - se, ymax= mean.delta.resid+se), width = 0.5, size=1.1) + geom_smooth(work.set3c, mapping = aes(x = lunar.day, y = delta.resid), stat="smooth",position="identity", method = "auto", level= 0.95, size=1.1, colour="black") + ylab("Residual growth (mm/d)\n") + xlab("Lunar day") + #relabelling x and y axes ggtitle(" (c) Larval ages: 3-10d")+ scale_x_continuous(limits = c(-1,30), breaks = seq(from = 0, to = 30, by = 5)) + theme_classic (base_size = 27) + theme(axis.title.y = element_text(margin = margin(t = 0, r = 28, b = 0, l = 0)), plot.title = element_text(colour="black", size=rel(1), hjust=0.05, vjust=-3)) figA4c work.set3d <- work.set %>% filter(!is.na(delta.resid)) %>% filter(larval.age>=10) %>% filter(larval.age<=20) work.set3d1 <- work.set3d %>% group_by(lunar.day) %>% summarise(mean.delta.resid = mean(delta.resid), n=n(), sd=sd(delta.resid), se=sd/sqrt(n)) figA4d <- ggplot(work.set3d1, aes(x = lunar.day, y = mean.delta.resid)) + geom_point(size = 5, alpha = .8) + geom_errorbar(aes(ymin = mean.delta.resid - se, ymax= mean.delta.resid+se), width = 0.5, size=1.1) + geom_smooth(work.set3d, mapping = aes(x = lunar.day, y = delta.resid), stat="smooth",position="identity", method = "auto", level= 0.95, size=1.1, colour="black") + ylab("Residual growth (mm/d)\n") + xlab("Lunar day") + ggtitle(" (d) Larval ages: 10-20d")+ scale_x_continuous(limits = c(-1,30), breaks = seq(from = 0, to = 30, by = 5)) + theme_classic (base_size = 27) + theme(axis.title.y = element_text(margin = margin(t = 0, r = 28, b = 0, l = 0)), plot.title = element_text(colour="black", size=rel(1), hjust=0.05, vjust=-3)) figA4d work.set3e <- work.set %>% filter(!is.na(delta.resid)) %>% filter(larval.age>=20) %>% filter(larval.age<=30) work.set3e1 <- work.set3e %>% group_by(lunar.day) %>% summarise(mean.delta.resid = mean(delta.resid), n=n(), sd=sd(delta.resid), se=sd/sqrt(n)) figA4e <- ggplot(work.set3e1, aes(x = lunar.day, y = mean.delta.resid)) + geom_point(size = 5, alpha = .8) + geom_errorbar(aes(ymin = mean.delta.resid - se, ymax= mean.delta.resid+se), width = 0.5, size=1.1) + geom_smooth(work.set3e, mapping = aes(x = lunar.day, y = delta.resid), stat="smooth",position="identity", method = "auto", level= 0.95, size=1.1, colour="black") + ylab("Residual growth (mm/d)\n") + xlab("Lunar day") + ggtitle(" (e) Larval ages: 20-30d")+ scale_x_continuous(limits = c(-1,30), breaks = seq(from = 0, to = 30, by = 5)) + theme_classic (base_size = 27) + theme(axis.title.y = element_text(margin = margin(t = 0, r = 28, b = 0, l = 0)), plot.title = element_text(colour="black", size=rel(1), hjust=0.05, vjust=-3)) figA4e work.set3f <- work.set %>% filter(!is.na(delta.resid)) %>% filter(larval.age>=30) %>% filter(larval.age<=40) work.set3f1 <- work.set3f %>% group_by(lunar.day) %>% summarise(mean.delta.resid = mean(delta.resid), n=n(), sd=sd(delta.resid), se=sd/sqrt(n)) figA4f <- ggplot(work.set3f1, aes(x = lunar.day, y = mean.delta.resid)) + geom_point(size = 5, alpha = .8) + geom_errorbar(aes(ymin = mean.delta.resid - se, ymax= mean.delta.resid+se), width = 0.5, size=1.1) + geom_smooth(work.set3f, mapping = aes(x = lunar.day, y = delta.resid), stat="smooth",position="identity", method = "auto", level= 0.95, size=1.1, colour="black") + ylab("Residual growth (mm/d)\n") + xlab("Lunar day") + ggtitle(" (f) Larval ages: 30-40d")+ scale_x_continuous(limits = c(-1,30), breaks = seq(from = 0, to = 30, by = 5)) + theme_classic (base_size = 27) + theme(axis.title.y = element_text(margin = margin(t = 0, r = 28, b = 0, l = 0)), plot.title = element_text(colour="black", size=rel(1), hjust=0.05, vjust=-3)) figA4f work.set3g <- work.set %>% filter(!is.na(delta.resid)) %>% filter(larval.age>=40) %>% filter(larval.age<=50) work.set3g1 <- work.set3g %>% group_by(lunar.day) %>% summarise(mean.delta.resid = mean(delta.resid), n=n(), sd=sd(delta.resid), se=sd/sqrt(n)) figA4g <- ggplot(work.set3g1, aes(x = lunar.day, y = mean.delta.resid)) + geom_point(size = 5, alpha = .8) + geom_errorbar(aes(ymin = mean.delta.resid - se, ymax= mean.delta.resid+se), width = 0.5, size=1.1) + geom_smooth(work.set3g, mapping = aes(x = lunar.day, y = delta.resid), stat="smooth",position="identity", method = "auto", level= 0.95, size=1.1, colour="black") + ylab("Residual growth (mm/d)\n") + xlab("Lunar day") + ggtitle(" (g) Larval ages: 40-50d")+ scale_x_continuous(limits = c(-1,30), breaks = seq(from = 0, to = 30, by = 5)) + theme_classic (base_size = 27) + theme(axis.title.y = element_text(margin = margin(t = 0, r = 28, b = 0, l = 0)), plot.title = element_text(colour="black", size=rel(1), hjust=0.05, vjust=-3)) figA4g work.set3h <- work.set %>% filter(!is.na(delta.resid)) %>% filter(larval.age>=50) %>% filter(larval.age<=60) work.set3h1 <- work.set3h %>% group_by(lunar.day) %>% summarise(mean.delta.resid = mean(delta.resid), n=n(), sd=sd(delta.resid), se=sd/sqrt(n)) figA4h <- ggplot(work.set3h1, aes(x = lunar.day, y = mean.delta.resid)) + geom_point(size = 5, alpha = .8) + geom_errorbar(aes(ymin = mean.delta.resid - se, ymax= mean.delta.resid+se), width = 0.5, size=1.1) + geom_smooth(work.set3h, mapping = aes(x = lunar.day, y = delta.resid), stat="smooth",position="identity", method = "auto", level= 0.95, size=1.1, colour="black") + ylab("Residual growth (mm/d)\n") + xlab("Lunar day") + ggtitle(" (h) Larval ages: 50-60d")+ scale_x_continuous(limits = c(-1,30), breaks = seq(from = 0, to = 30, by = 5)) + theme_classic (base_size = 27) + theme(axis.title.y = element_text(margin = margin(t = 0, r = 28, b = 0, l = 0)), plot.title = element_text(colour="black", size=rel(1), hjust=0.05, vjust=-3)) figA4h figA4<- grid.arrange(figA4a, figA4b, figA4c, figA4d, figA4e, figA4f, figA4g, figA4h, ncol = 2, nrow = 4) # ggsave(plot = figA4, width=20, height=32, dpi=300,filename = "figures\\Fig A3.pdf") # ggsave(plot = figA4, width=20, height=32, dpi=300,filename = "figures\\Fig A3.png") ####################################################### ####################################################### # (6) Explore potential explanations for lunar periodicity in residual growth using a secondary analysis #importing file that contains data describing the proportion of 18:00-0:00 when moon is above horizon ("proportion.first.half"), and proportion of 0:00-06:00 when moon is above horizon ("proportion.second.half"), excel file of same name shows how this file was constructed moonlight2<-read.csv("moonlight_ancillary_fig_3.csv") moonlight2 <- moonlight2 %>% select(date, proportion.first.half, proportion.second.half)%>% mutate(date=dmy(date)) #correcting date format work.set2.sources <- left_join(work.set2,moonlight2,by ="date" ) %>% mutate(theta.bd=(360*LD.birth)/28.5, #Adding parameters to facilitate periodic regression; 28.5d because numbering starts at 0 sin.theta.bd=sin(theta.bd*pi/180), cos.theta.bd=cos(theta.bd*pi/180), logit.cloud.3x3=logit(cloud.cover.airs.nighttime.3x3), moonlight.early=proportion.first.half*illumination, moonlight.late=proportion.second.half*illumination) glimpse(work.set2.sources) ggplot(data = work.set2.sources, aes(x = logit.cloud.3x3 )) + geom_histogram(binwidth=.025, colour="black") #logit transform normalizes data #################### # The model (with) cloud cover) # Because residuals of a base model exhibit patterns of autocorrelation, will explore AR1 and ARMA error structure # lmer (in package lme4) does not easily allow for modeling spatial/temporal autocorrelation ##################################### m.sources.lme.ar1<- lme(delta.resid ~ moonlight.early + moonlight.early/logit.cloud.3x3 + moonlight.late + moonlight.late/logit.cloud.3x3 , random = ~1 | fish.id, work.set2.sources, correlation=corAR1(form = ~date),method="REML") m.sources.lme.arma<- lme(delta.resid ~ moonlight.early + moonlight.early/logit.cloud.3x3 + moonlight.late + moonlight.late/logit.cloud.3x3, random = ~1 | fish.id, work.set2.sources, correlation=corARMA(form = ~date, p = 1, q = 1),method="REML") AICc(m.sources.lme.ar1,m.sources.lme.arma) #ARMA model is better plot(m.sources.lme.ar1) qqnorm(resid(m.sources.lme.ar1)) qqline(resid(m.sources.lme.ar1)) acf(residuals(m.sources.lme.ar1, type = "normalized", lag=30)) # autocorrelation in residuals has been partially eliminated plot(m.sources.lme.arma) qqnorm(resid(m.sources.lme.arma)) qqline(resid(m.sources.lme.arma)) acf(residuals(m.sources.lme.arma, type = "normalized", lag=30)) # autocorrelation in residuals has been eliminated anova(m.sources.lme.arma) summary(m.sources.lme.arma) ##################################### #################### #Figure 3 (plots of predicted values) e1.m.sources <- predictorEffect("logit.cloud.3x3", m.sources.lme.arma, xlevels=list(illumination=c(0,0.5,1), proportion.first.half=c(0,1), proportion.second.half=c(0,1))) #using effects package to generate a data set of predicted values for "cloud.fraction" conditional on a set of values for other x variables e1.m.sources.df <- data.frame(e1.m.sources) e1.m.sources.df <- e1.m.sources.df %>% mutate(cloud.fraction=(exp(logit.cloud.3x3+0.025)/(1+exp(logit.cloud.3x3+0.025)))) new.moon <- e1.m.sources.df%>% #subsetting the dataframe to correspond to x values for particular moon phases filter(moonlight.early==0, moonlight.late==0) first.quarter <- e1.m.sources.df%>% filter(moonlight.early==0.5,moonlight.late==0) full.moon <- e1.m.sources.df%>% filter(moonlight.early==1,moonlight.late==1) last.quarter <- e1.m.sources.df%>% filter(moonlight.early==0, moonlight.late==0.5) moon.x1 <- -0.17 moon.x2 <- 0.105 moon.x3 <- 0.38 moon.y <- -0.25 moon.scale<-0.17 cloud.x2 <- 0.06 cloud.x3.1 <- 0.34 cloud.x3.2 <- 0.41 cloud.y1 <- -0.27 cloud.y2 <- -0.255 cloud.scale<-0.1 cloud.scale2<-0.117 fig3a<-ggplot(new.moon, aes(x= cloud.fraction, y = fit))+ #geom_smooth()+ geom_smooth(aes(ymin = lower, ymax = upper), stat='identity', size=1.1,colour="black")+ geom_smooth(aes(y=0), stat='identity', size=0.5, colour="black", linetype="dashed") + ggtitle(" (a) New moon")+ ylab("Residual growth (mm/d)\n") + xlab("Cloud cover") + scale_x_continuous(limits = c(0,1), breaks = seq(from = 0, to = 1, by = 0.2)) + scale_y_continuous(limits = c(-0.03,0.02), breaks = seq(from = -0.03, to = 0.02, by = 0.01)) + theme_classic (base_size = 34) + theme(axis.title.y = element_text(margin = margin(t = 0, r = -30, b = 0, l = -10)), plot.title = element_text(colour="black", size=rel(1), hjust=0.05, vjust=-3), axis.title.x = element_text(margin = margin(t = 0, r = 0, b = 0, l = 0))) fig3a <- ggdraw(fig3a) + draw_image("figures\\moon_icon\\new_moon.PNG", x= moon.x1, y = moon.y, scale=moon.scale) + draw_image("figures\\moon_icon\\new_moon.PNG", x= moon.x2, y = moon.y, scale=moon.scale) + draw_image("figures\\moon_icon\\new_moon.PNG", x= moon.x3, y = moon.y, scale=moon.scale) + draw_image("figures\\moon_icon\\cloud-b.png", x= cloud.x2, y = cloud.y1, scale=cloud.scale2) + draw_image("figures\\moon_icon\\cloud-b.png", x= cloud.x3.1, y = cloud.y1, scale=cloud.scale2) + draw_image("figures\\moon_icon\\cloud-b.png", x= cloud.x3.2, y = cloud.y2, scale=cloud.scale2) + draw_image("figures\\moon_icon\\cloud.png", x= cloud.x2, y = cloud.y1, scale=cloud.scale) + draw_image("figures\\moon_icon\\cloud.png", x= cloud.x3.1, y = cloud.y1, scale=cloud.scale) + draw_image("figures\\moon_icon\\cloud.png", x= cloud.x3.2, y = cloud.y2, scale=cloud.scale) fig3a fig3b<-ggplot(first.quarter, aes(x= cloud.fraction, y = fit))+ #geom_smooth()+ geom_smooth(aes(ymin = lower, ymax = upper), stat='identity', size=1.1,colour="black")+ geom_smooth(aes(y=0), stat='identity', size=0.5, colour="black", linetype="dashed") + ggtitle(" (b) First quarter moon")+ ylab("Residual growth (mm)\n") + xlab("Cloud cover") + scale_x_continuous(limits = c(0,1), breaks = seq(from = 0, to = 1, by = 0.2)) + scale_y_continuous(limits = c(-0.03,0.02), breaks = seq(from = -0.03, to = 0.02, by = 0.01)) + theme_classic (base_size = 34) + theme(axis.title.y = element_text(colour="white", margin = margin(t = 0, r = -30, b = 0, l = -10)), plot.title = element_text(colour="black", size=rel(1), hjust=0.05, vjust=-3), axis.title.x = element_text(margin = margin(t = 0, r = 0, b = 0, l = 0))) fig3b <- ggdraw(fig3b) + draw_image("figures\\moon_icon\\first_quarter.PNG", x= moon.x1, y = moon.y, scale=moon.scale) + draw_image("figures\\moon_icon\\first_quarter.PNG", x= moon.x2, y = moon.y, scale=moon.scale) + draw_image("figures\\moon_icon\\first_quarter.PNG", x= moon.x3, y = moon.y, scale=moon.scale) + draw_image("figures\\moon_icon\\cloud-b.png", x= cloud.x2, y = cloud.y1, scale=cloud.scale2) + draw_image("figures\\moon_icon\\cloud-b.png", x= cloud.x3.1, y = cloud.y1, scale=cloud.scale2) + draw_image("figures\\moon_icon\\cloud-b.png", x= cloud.x3.2, y = cloud.y2, scale=cloud.scale2) + draw_image("figures\\moon_icon\\cloud.png", x= cloud.x2, y = cloud.y1, scale=cloud.scale) + draw_image("figures\\moon_icon\\cloud.png", x= cloud.x3.1, y = cloud.y1, scale=cloud.scale) + draw_image("figures\\moon_icon\\cloud.png", x= cloud.x3.2, y = cloud.y2, scale=cloud.scale) fig3b fig3c<-ggplot(full.moon, aes(x= cloud.fraction, y = fit))+ #geom_smooth()+ geom_smooth(aes(ymin = lower, ymax = upper), stat='identity', size=1.1,colour="black")+ geom_smooth(aes(y=0), stat='identity', size=0.5, colour="black", linetype="dashed") + ggtitle(" (c) Full moon")+ ylab("Residual growth (mm)\n") + xlab("Cloud cover") + scale_x_continuous(limits = c(0,1), breaks = seq(from = 0, to = 1, by = 0.2)) + scale_y_continuous(limits = c(-0.03,0.02), breaks = seq(from = -0.03, to = 0.02, by = 0.01)) + theme_classic (base_size = 34) + theme(axis.title.y = element_text(colour="white", margin = margin(t = 0, r = -30, b = 0, l = -10)), plot.title = element_text(colour="black", size=rel(1), hjust=0.05, vjust=-3), axis.title.x = element_text(margin = margin(t = 0, r = 0, b = 0, l = 0))) fig3c <- ggdraw(fig3c) + draw_image("figures\\moon_icon\\full_moon.PNG", x= moon.x1, y = moon.y, scale=moon.scale) + draw_image("figures\\moon_icon\\full_moon.PNG", x= moon.x2, y = moon.y, scale=moon.scale) + draw_image("figures\\moon_icon\\full_moon.PNG", x= moon.x3, y = moon.y, scale=moon.scale) + draw_image("figures\\moon_icon\\cloud-b.png", x= cloud.x2, y = cloud.y1, scale=cloud.scale2) + draw_image("figures\\moon_icon\\cloud-b.png", x= cloud.x3.1, y = cloud.y1, scale=cloud.scale2) + draw_image("figures\\moon_icon\\cloud-b.png", x= cloud.x3.2, y = cloud.y2, scale=cloud.scale2) + draw_image("figures\\moon_icon\\cloud.png", x= cloud.x2, y = cloud.y1, scale=cloud.scale) + draw_image("figures\\moon_icon\\cloud.png", x= cloud.x3.1, y = cloud.y1, scale=cloud.scale) + draw_image("figures\\moon_icon\\cloud.png", x= cloud.x3.2, y = cloud.y2, scale=cloud.scale) fig3c fig3d<-ggplot(last.quarter, aes(x= cloud.fraction, y = fit))+ #geom_smooth()+ geom_smooth(aes(ymin = lower, ymax = upper), stat='identity', size=1.1,colour="black")+ geom_smooth(aes(y=0), stat='identity', size=0.5, colour="black", linetype="dashed") + ggtitle(" (d) Last quarter moon")+ ylab("Residual growth (mm)\n") + xlab("Cloud cover") + scale_x_continuous(limits = c(0,1), breaks = seq(from = 0, to = 1, by = 0.2)) + scale_y_continuous(limits = c(-0.03,0.02), breaks = seq(from = -0.03, to = 0.02, by = 0.01)) + theme_classic (base_size = 34) + theme(axis.title.y = element_text(colour="white", margin = margin(t = 0, r = -30, b = 0, l = -10)), plot.title = element_text(colour="black", size=rel(1), hjust=0.05, vjust=-3), axis.title.x = element_text(margin = margin(t = 0, r = 0, b = 0, l = 0))) fig3d <- ggdraw(fig3d) + draw_image("figures\\moon_icon\\third_quarter.PNG", x= moon.x1, y = moon.y, scale=moon.scale) + draw_image("figures\\moon_icon\\third_quarter.PNG", x= moon.x2, y = moon.y, scale=moon.scale) + draw_image("figures\\moon_icon\\third_quarter.PNG", x= moon.x3, y = moon.y, scale=moon.scale) + draw_image("figures\\moon_icon\\cloud-b.png", x= cloud.x2, y = cloud.y1, scale=cloud.scale2) + draw_image("figures\\moon_icon\\cloud-b.png", x= cloud.x3.1, y = cloud.y1, scale=cloud.scale2) + draw_image("figures\\moon_icon\\cloud-b.png", x= cloud.x3.2, y = cloud.y2, scale=cloud.scale2) + draw_image("figures\\moon_icon\\cloud.png", x= cloud.x2, y = cloud.y1, scale=cloud.scale) + draw_image("figures\\moon_icon\\cloud.png", x= cloud.x3.1, y = cloud.y1, scale=cloud.scale) + draw_image("figures\\moon_icon\\cloud.png", x= cloud.x3.2, y = cloud.y2, scale=cloud.scale) fig3d fig3<- grid.arrange(fig3a, fig3b, fig3c, fig3d, ncol = 4, nrow = 1) #ggsave(plot = fig3, width=30, height=7, dpi=300,filename = "figures\\Fig 3.pdf") #ggsave(plot = fig3, width=30, height=7, dpi=300,filename = "figures\\Fig 3.png") #END #######################################################################