library(ggplot2) library(ggpubr) library(mgcv) library(visreg) library(strucchange) library(cmocean) library(segmented) library(DLMtool) #-----------------------------------------------------# # 1. Merging bull shark isotope data with metadata #### #-----------------------------------------------------# df <- read.csv("Input/BS_isotopes.csv") df <- df[-101,] # Remove outlier! df$Name <- as.character(df$Name) df.bio <- read.csv("Input/SIA_samples.csv") df.bio$ID <- paste0(df.bio$SampleLabel, "-", substr(df.bio$Sample, 1,1)) df.bio$Sex <- as.character(df.bio$Sex) df.bio$Sample <- as.character(df.bio$Sample) df$Sex <- NA df$TL <- NA df$Weight.kg. <- NA df$Sample <- NA for (i in 1:nrow(df)) { aux <- df.bio[which(df.bio$ID == df$Name[i]), ] df$Sex[i] <- aux$Sex[1] df$TL[i] <- aux$TL.cm[1] df$Weight.kg.[i] <- aux$Weight[1] df$Sample[i] <- aux$Sample[1] } # Remove NA missing from 103 cm Male Liver df <- df[-which(is.na(df$Weight.kg.) == TRUE), ] # Add relative weight variable: df$Weight.rel <- (((df$Weight.kg. * 100) / 2.54) - 100) / 100 # Age: max = 5 years t <- seq(1, 6.5, 0.5) # First year size t.save <- NULL for (i in 1:length(t)) { aux <- 256.4 * (1 - exp(-0.139 * (t[i] + 1.9))) # Cruz-Martinez et al 2005 (both sexes combined) t.save <- c(t.save, aux) } # Inverted equation to find out ages from size: df$Age.n <- NA for (i in 1:nrow(df)) { df$Age.n[i] <- iVB(t0 = -1.9, K = 0.139, Linf = 256.4, L = df$TL[i]) } df$Age <- "Juvenile" df$Age[df$TL <= t.save[1]] <- "YOY" # Plot relationships mod.tl <- gam(Weight.kg. ~ s(TL, k = 4), family = gaussian, data = df) mod.tl2 <- lm(Weight.kg. ~ TL, data = df) mod.tl2 <- segmented(mod.tl2, seg.Z = ~TL) fit <- numeric(length(df$TL)) * NA fit[complete.cases(rowSums(cbind(df$Weight.kg., df$TL)))] <- broken.line(mod.tl2)$fit data1 <- data.frame(W = df$Weight.kg., TL = df$TL, fit = fit) visreg(mod.tl, partial = FALSE, gg = TRUE, rug = FALSE, type = "conditional", line = c(lwd = 2, col = "gray"), fill.par = c(fill = NA), xlab = "Total length (cm)", ylab = "Weight (Kg)") + theme_bw() + geom_line(data = data1, aes(x = TL, y = fit), color = 'blue') + geom_point(data = df, aes(x = TL, y = Weight.kg., color = Age, shape = Sex)) + scale_color_manual(values = c("black", "red")) + scale_x_continuous(sec.axis = sec_axis(~ . * 1, breaks = c(75, t.save), labels = c("YOY", "1", "", "2", "", "3", "", "4", "", "5", "", "6", ""), name = "Age (years)")) + scale_y_continuous(sec.axis = sec_axis(~ . * 1, name = "Relative Increase in Weight (RIW)", breaks = c(5.10, 14.60, 24.50, 35.40), labels = c("1.00", "4.75", "8.65", "13.00"))) + geom_vline(xintercept = t.save[1], lty = "dashed", col = "red") + theme(legend.position = "top") + labs(shape = "", color = "") + annotate("text", x = 110, y = 36, label = expression(paste("TL = 256.4(1-e"^"-0.139(Age+1.9)", ")")), size = 4) ggsave(filename = "Output/Length_weight.png", width = 18, height = 14, units = "cm") #------------------------------------# # 2. Exploratory analysis of data #### #------------------------------------# df.m <- subset(df, Sample == "Muscle") df.l <- subset(df, Sample == "Liver") # N/C/S x TL (Muscle) plot1 <- ggplot() + geom_point(data = df.m, aes(x = TL, y = d15NAIR, color = Sex)) + theme_classic() + labs(x = "TL (cm)", y = "N", title = "Muscle") plot2 <- ggplot() + geom_point(data = df.m, aes(x = TL, y = d13CVPDB, color = Sex)) + theme_classic() + labs(x = "TL (cm)", y = "C", title = "Muscle") plot3 <- ggplot() + geom_point(data = df.m, aes(x = TL, y = d34SVCDT, color = Sex)) + theme_classic() + labs(x = "TL (cm)", y = "S", title = "Muscle") ggarrange(plot1, plot2, plot3, ncol = 3) ggsave(filename = "Output/N_C_TL_muscle.jpeg", width = 30, height = 10, units = "cm") # N/C x TL (Liver) plot1 <- ggplot() + geom_point(data = df.l, aes(x = TL, y = d15NAIR, color = Sex)) + theme_classic() + labs(x = "TL (cm)", y = "N", title = "Liver") plot2 <- ggplot() + geom_point(data = df.l, aes(x = TL, y = d13CVPDB, color = Sex)) + theme_classic() + labs(x = "TL (cm)", y = "C", title = "Liver") plot3 <- ggplot() + geom_point(data = df.l, aes(x = TL, y = d34SVCDT, color = Sex)) + theme_classic() + labs(x = "TL (cm)", y = "S", title = "Liver") ggarrange(plot1, plot2, plot3, ncol = 3) ggsave(filename = "Output/N_C_TL_liver.jpeg", width = 30, height = 10, units = "cm") #----------------------# # 3. Analysis Start #### #----------------------# ### Comparing model types: # Muscle: # N: mod1.full <- gam(d15NAIR ~ s(TL, k = 4), family = gaussian, data = df.m) # GAM mod2.full <- gam(d15NAIR ~ s(Weight.rel, k = 4), family = gaussian, data = df.m) # GAM mod1.yoy <- gam(d15NAIR ~ s(TL, k = 4), subset = df.m$TL > t.save[1], family = gaussian, data = df.m) # GAM mod2.yoy <- gam(d15NAIR ~ s(Weight.rel, k = 4), subset = df.m$TL > t.save[1],family = gaussian, data = df.m) # GAM summary(mod1.full) summary(mod1.yoy) summary(mod2.full) summary(mod2.yoy) # C: mod1.full <- gam(d13CVPDB ~ s(TL, k = 4), family = gaussian, data = df.m) # GAM mod2.full <- gam(d13CVPDB ~ s(Weight.rel, k = 4), family = gaussian, data = df.m) # GAM mod1.yoy <- gam(d13CVPDB ~ s(TL, k = 4), subset = df.m$TL > t.save[1], family = gaussian, data = df.m) # GAM mod2.yoy <- gam(d13CVPDB ~ s(Weight.rel, k = 4), subset = df.m$TL > t.save[1],family = gaussian, data = df.m) # GAM summary(mod1.full) summary(mod1.yoy) summary(mod2.full) summary(mod2.yoy) # S: mod1.full <- gam(d34SVCDT ~ s(TL, k = 4), family = gaussian, data = df.m) # GAM mod2.full <- gam(d34SVCDT ~ s(Weight.rel, k = 4), family = gaussian, data = df.m) # GAM mod1.yoy <- gam(d34SVCDT ~ s(TL, k = 4), subset = df.m$TL > t.save[1], family = gaussian, data = df.m) # GAM mod2.yoy <- gam(d34SVCDT ~ s(Weight.rel, k = 4), subset = df.m$TL > t.save[1],family = gaussian, data = df.m) # GAM summary(mod1.full) summary(mod1.yoy) summary(mod2.full) summary(mod2.yoy) # Liver: # N: mod1.full <- gam(d15NAIR ~ s(TL, k = 4), family = gaussian, data = df.l) # GAM mod2.full <- gam(d15NAIR ~ s(Weight.rel, k = 4), family = gaussian, data = df.l) # GAM mod1.yoy <- gam(d15NAIR ~ s(TL, k = 4), subset = df.l$TL > t.save[1], family = gaussian, data = df.l) # GAM mod2.yoy <- gam(d15NAIR ~ s(Weight.rel, k = 4), subset = df.l$TL > t.save[1],family = gaussian, data = df.l) # GAM summary(mod1.full) summary(mod1.yoy) summary(mod2.full) summary(mod2.yoy) # C: mod1.full <- gam(d13CVPDB ~ s(TL, k = 4), family = gaussian, data = df.l) # GAM mod2.full <- gam(d13CVPDB ~ s(Weight.rel, k = 4), family = gaussian, data = df.l) # GAM mod1.yoy <- gam(d13CVPDB ~ s(TL, k = 4), subset = df.l$TL > t.save[1], family = gaussian, data = df.l) # GAM mod2.yoy <- gam(d13CVPDB ~ s(Weight.rel, k = 4), subset = df.l$TL > t.save[1],family = gaussian, data = df.l) # GAM summary(mod1.full) summary(mod1.yoy) summary(mod2.full) summary(mod2.yoy) # S: mod1.full <- gam(d34SVCDT ~ s(TL, k = 4), family = gaussian, data = df.l) # GAM mod2.full <- gam(d34SVCDT ~ s(Weight.rel, k = 4), family = gaussian, data = df.l) # GAM mod1.yoy <- gam(d34SVCDT ~ s(TL, k = 4), subset = df.l$TL > t.save[1], family = gaussian, data = df.l) # GAM mod2.yoy <- gam(d34SVCDT ~ s(Weight.rel, k = 4), subset = df.l$TL > t.save[1],family = gaussian, data = df.l) # GAM summary(mod1.full) summary(mod1.yoy) summary(mod2.full) summary(mod2.yoy) ### Ploting final models ## Muscle: # Nitrogren mod.gam.n.1 <- gam(d15NAIR ~ s(TL, k = 4), family = gaussian, data = df.m) # GAM mod.gam.n.2 <- gam(d15NAIR ~ s(Weight.rel, k = 4), family = gaussian, data = df.m) # GAM summary(mod.gam.n.1) summary(mod.gam.n.2) plot1 <- visreg(mod.gam.n.1, partial = FALSE, gg = TRUE, rug = TRUE, type = "conditional", line = c(lwd = 2, col = cmocean("haline")(4)[1]), fill.par = c(fill = rgb(col2rgb(cmocean('haline')(4)[1])[1,], col2rgb(cmocean('haline')(4)[1])[2,], col2rgb(cmocean('haline')(4)[1])[3,], max = 255, alpha = 125)), ylab = expression(paste(delta^15, "N")), ylim = c(5, 20), xlim = c(70, 180)) + theme_bw() + geom_point(data = df.m, aes(x = TL, y = d15NAIR, shape = Sex), col = cmocean('haline')(4)[1]) + theme(legend.position = "none") + scale_x_continuous(breaks = seq(75, 175, 25), name = "TL (cm)") plot2 <- visreg(mod.gam.n.2, partial = FALSE, gg = TRUE, rug = TRUE, type = "conditional", line = c(lwd = 2, col = cmocean("haline")(4)[1]), fill.par = c(fill = rgb(col2rgb(cmocean('haline')(4)[1])[1,], col2rgb(cmocean('haline')(4)[1])[2,], col2rgb(cmocean('haline')(4)[1])[3,], max = 255, alpha = 125)), xlab = "RIW", ylab = expression(paste(delta^15, "N")), ylim = c(5, 20), xlim = c(0, 15)) + theme_bw() + geom_point(data = df.m, aes(x = Weight.rel, y = d15NAIR, shape = Sex), col = cmocean('haline')(4)[1]) + theme(legend.position = "none") # Carbon mod.gam.c.1 <- gam(d13CVPDB ~ s(TL, k = 4), family = gaussian, data = df.m) # GAM mod.gam.c.2 <- gam(d13CVPDB ~ s(Weight.rel, k = 4), family = gaussian, data = df.m) # GAM summary(mod.gam.c.1) summary(mod.gam.c.2) plot3 <- visreg(mod.gam.c.1, partial = FALSE, gg = TRUE, rug = TRUE, type = "conditional", line = c(lwd = 2, col = cmocean("haline")(4)[2]), fill.par = c(fill = rgb(col2rgb(cmocean('haline')(4)[2])[1,], col2rgb(cmocean('haline')(4)[2])[2,], col2rgb(cmocean('haline')(4)[2])[3,], max = 255, alpha = 125)), xlab = "Total length (cm)", ylab = expression(paste(delta^13, "C")), ylim = c(5, 20), xlim = c(70, 180)) + theme_bw() + geom_point(data = df.m, aes(x = TL, y = d13CVPDB, shape = Sex), col = cmocean('haline')(4)[2]) + theme(legend.position = "none") + scale_x_continuous(breaks = seq(75, 175, 25), name = "TL (cm)") plot4 <- visreg(mod.gam.c.2, partial = FALSE, gg = TRUE, rug = TRUE, type = "conditional", line = c(lwd = 2, col = cmocean("haline")(4)[2]), fill.par = c(fill = rgb(col2rgb(cmocean('haline')(4)[2])[1,], col2rgb(cmocean('haline')(4)[2])[2,], col2rgb(cmocean('haline')(4)[2])[3,], max = 255, alpha = 125)), xlab = "RIW", ylab = expression(paste(delta^13, "C")), ylim = c(5, 20), xlim = c(0, 15)) + theme_bw() + geom_point(data = df.m, aes(x = Weight.rel, y = d13CVPDB, shape = Sex), col = cmocean('haline')(4)[2]) + theme(legend.position = "none") # Sulfur mod.gam.s.1 <- gam(d34SVCDT ~ s(TL, k = 4), family = gaussian, data = df.m) # GAM mod.gam.s.2 <- gam(d34SVCDT ~ s(Weight.rel, k = 4), family = gaussian, data = df.m) # GAM summary(mod.gam.s.1) summary(mod.gam.s.2) plot5 <- visreg(mod.gam.s.1, partial = FALSE, gg = TRUE, rug = TRUE, type = "conditional", line = c(lwd = 2, col = cmocean("haline")(4)[3]), fill.par = c(fill = rgb(col2rgb(cmocean('haline')(4)[3])[1,], col2rgb(cmocean('haline')(4)[3])[2,], col2rgb(cmocean('haline')(4)[3])[3,], max = 255, alpha = 125)), xlab = "Total length (cm)", ylab = expression(paste(delta^34, "S")), ylim = c(5, 20), xlim = c(70, 180)) + theme_bw() + geom_point(data = df.m, aes(x = TL, y = d34SVCDT, shape = Sex), col = cmocean('haline')(4)[3]) + theme(legend.position = "none") + scale_x_continuous(breaks = seq(75, 175, 25), name = "TL (cm)") plot6 <- visreg(mod.gam.s.2, partial = FALSE, gg = TRUE, rug = TRUE, type = "conditional", line = c(lwd = 2, col = cmocean("haline")(4)[3]), fill.par = c(fill = rgb(col2rgb(cmocean('haline')(4)[3])[1,], col2rgb(cmocean('haline')(4)[3])[2,], col2rgb(cmocean('haline')(4)[3])[3,], max = 255, alpha = 125)), xlab = "RIW", ylab = expression(paste(delta^34, "S")), ylim = c(5, 20), xlim = c(0, 15)) + theme_bw() + geom_point(data = df.m, aes(x = Weight.rel, y = d34SVCDT, shape = Sex), col = cmocean('haline')(4)[3]) + theme(legend.position = "none") ## Liver: # Nitrogren mod.gam.n.1 <- gam(d15NAIR ~ s(TL, k = 4), subset = TL > 80, family = gaussian, data = df.l) # GAM mod.gam.n.2 <- gam(d15NAIR ~ s(Weight.rel, k = 4), subset = TL > 80, family = gaussian, data = df.l) # GAM summary(mod.gam.n.1) summary(mod.gam.n.2) plot7 <- visreg(mod.gam.n.1, partial = FALSE, gg = TRUE, rug = TRUE, type = "conditional", line = c(lwd = 2, col = cmocean("haline")(4)[1]), fill.par = c(fill = rgb(col2rgb(cmocean('haline')(4)[1])[1,], col2rgb(cmocean('haline')(4)[1])[2,], col2rgb(cmocean('haline')(4)[1])[3,], max = 255, alpha = 125)), xlab = "Total length (cm)", ylab = expression(paste(delta^15, "N")), ylim = c(5, 20), xlim = c(70, 180)) + theme_bw() + geom_point(data = df.l, aes(x = TL, y = d15NAIR, shape = Sex), col = cmocean('haline')(4)[1]) + theme(legend.position = "none") + scale_x_continuous(breaks = seq(75, 175, 25), name = "TL (cm)") plot8 <- visreg(mod.gam.n.2, partial = FALSE, gg = TRUE, rug = TRUE, type = "conditional", line = c(lwd = 2, col = cmocean("haline")(4)[1]), fill.par = c(fill = rgb(col2rgb(cmocean('haline')(4)[1])[1,], col2rgb(cmocean('haline')(4)[1])[2,], col2rgb(cmocean('haline')(4)[1])[3,], max = 255, alpha = 125)), xlab = "RIW", ylab = expression(paste(delta^15, "N")), ylim = c(5, 20), xlim = c(0, 15)) + theme_bw() + geom_point(data = df.l, aes(x = Weight.rel, y = d15NAIR, shape = Sex), col = cmocean('haline')(4)[1]) + theme(legend.position = "none") # Carbon mod.gam.c.1 <- gam(d13CVPDB ~ s(TL, k = 4), subset = TL > 80, family = gaussian, data = df.l) # GAM mod.gam.c.2 <- gam(d13CVPDB ~ s(Weight.rel, k = 4), subset = TL > 80, family = gaussian, data = df.l) # GAM summary(mod.gam.c.1) summary(mod.gam.c.2) plot9 <- visreg(mod.gam.c.1, partial = FALSE, gg = TRUE, rug = TRUE, type = "conditional", line = c(lwd = 2, col = cmocean("haline")(4)[2]), fill.par = c(fill = rgb(col2rgb(cmocean('haline')(4)[2])[1,], col2rgb(cmocean('haline')(4)[2])[2,], col2rgb(cmocean('haline')(4)[2])[3,], max = 255, alpha = 125)), xlab = "Total length (cm)", ylab = expression(paste(delta^13, "C")), ylim = c(5, 20), xlim = c(70, 180)) + theme_bw() + geom_point(data = df.l, aes(x = TL, y = d13CVPDB, shape = Sex), col = cmocean('haline')(4)[2]) + theme(legend.position = "none") + scale_x_continuous(breaks = seq(75, 175, 25), name = "TL (cm)") plot10 <- visreg(mod.gam.c.2, partial = FALSE, gg = TRUE, rug = TRUE, type = "conditional", line = c(lwd = 2, col = cmocean("haline")(4)[2]), fill.par = c(fill = rgb(col2rgb(cmocean('haline')(4)[2])[1,], col2rgb(cmocean('haline')(4)[2])[2,], col2rgb(cmocean('haline')(4)[2])[3,], max = 255, alpha = 125)), xlab = "RIW", ylab = expression(paste(delta^13, "C")), ylim = c(5, 20), xlim = c(0, 15)) + theme_bw() + geom_point(data = df.l, aes(x = Weight.rel, y = d13CVPDB, shape = Sex), col = cmocean('haline')(4)[2]) + theme(legend.position = "none") # Sulfur mod.gam.s.1 <- gam(d34SVCDT ~ s(TL, k = 4), subset = TL > 80, family = gaussian, data = df.l) # GAM mod.gam.s.2 <- gam(d34SVCDT ~ s(Weight.rel, k = 4), subset = TL > 80, family = gaussian, data = df.l) # GAM summary(mod.gam.s.1) summary(mod.gam.s.2) plot11 <- visreg(mod.gam.s.1, partial = FALSE, gg = TRUE, rug = TRUE, type = "conditional", line = c(lwd = 2, col = cmocean("haline")(4)[3]), fill.par = c(fill = rgb(col2rgb(cmocean('haline')(4)[3])[1,], col2rgb(cmocean('haline')(4)[3])[2,], col2rgb(cmocean('haline')(4)[3])[3,], max = 255, alpha = 125)), xlab = "Total length (cm)", ylab = expression(paste(delta^34, "S")), ylim = c(5, 20), xlim = c(70, 180)) + theme_bw() + geom_point(data = df.l, aes(x = TL, y = d34SVCDT, shape = Sex), col = cmocean('haline')(4)[3]) + theme(legend.position = "none") + scale_x_continuous(breaks = seq(75, 175, 25), name = "TL (cm)") plot12 <- visreg(mod.gam.s.2, partial = FALSE, gg = TRUE, rug = TRUE, type = "conditional", line = c(lwd = 2, col = cmocean("haline")(4)[3]), fill.par = c(fill = rgb(col2rgb(cmocean('haline')(4)[3])[1,], col2rgb(cmocean('haline')(4)[3])[2,], col2rgb(cmocean('haline')(4)[3])[3,], max = 255, alpha = 125)), xlab = "RIW", ylab = expression(paste(delta^34, "S")), ylim = c(5, 20), xlim = c(0, 15)) + theme_bw() + geom_point(data = df.l, aes(x = Weight.rel, y = d34SVCDT, shape = Sex), col = cmocean('haline')(4)[3]) + theme(legend.position = "none") # PLOT FINAL MODELS ggarrange(plot1, plot2, plot7, plot8, plot3, plot4, plot9, plot10, plot5, plot6, plot11, plot12, ncol = 4, nrow = 3, labels = c("A", "", "B", "", "", "", "", "", "", "", "", "")) ggsave("Output/Plot_reg.png", width = 25, height = 17, units = "cm") ### Find size breakpoints: segmented regressions ## Muscle # Nitrogen mod.base <- lm(d15NAIR ~ TL, data = df.m) mod.n <- segmented(mod.base, seg.Z = ~TL) # Carbon mod.base <- lm(d13CVPDB ~ TL, data = df.m) mod.c <- segmented(mod.base, seg.Z = ~TL) # Sulphur mod.base <- lm(d34SVCDT ~ TL, data = df.m) mod.s <- segmented(mod.base, seg.Z = ~TL) ## Liver # Nitrogen mod.base <- lm(d15NAIR ~ TL, data = df.l, subset = TL > 80) mod.n2 <- segmented(mod.base, seg.Z = ~TL) # Carbon mod.base <- lm(d13CVPDB ~ TL, data = df.l, subset = TL > 80) mod.c2 <- segmented(mod.base, seg.Z = ~TL) # Sulphur mod.base <- lm(d34SVCDT ~ TL, data = df.l, subset = TL > 80) mod.s2 <- segmented(mod.base, seg.Z = ~TL) ## Plot segmented regressions jpeg("Output/Break_muscle.png", width=22, height=22, units="cm", pointsize=12, quality=500,bg="white", res=500) par(mfrow=c(2, 1)) # Muscle par(mar = c(4, 12, 4, 1)) plot(df.m$d15NAIR ~ df.m$TL, type = "n", xlab = "", ylab = "", axes = FALSE, xlim = c(70, 175)) polygon(x = c(-100, -100, 135.3, 135.3), y = c(-1000, 1000, 1000, -10000), col = "gray90", border = NA) par(new = TRUE) plot(df.m$d15NAIR ~ df.m$TL, pch = 15, cex = 0.8, xlab = "", ylab = "", axes = FALSE, col = cmocean('haline')(4)[1], xlim = c(70, 175)) axis(1, at = seq(70, 170, 10)) axis(2, las = 1) axis(3, at = c(70, t.save), labels = c("0","1", "", "2", "", "3", "", "4", "", "5", "", "6", "")) mtext(text = "Age (years)", side = 3, line = 2.5) box() mtext(text = expression(paste(delta^15, "N")), side = 2, line = 2.5) par(new = TRUE) plot(df.m$d13CVPDB ~ df.m$TL, pch = 16, cex = 0.8, xlab = "", ylab = "", axes = FALSE, col = cmocean('haline')(4)[2], xlim = c(70, 175)) plot(mod.c, add=TRUE, link=FALSE, lwd=2, col = cmocean('haline')(4)[2]) axis(2, line = 4, las = 1) mtext(text = expression(paste(delta^13, "C")), side = 2, line = 6.5) par(new = TRUE) plot(df.m$d34SVCDT ~ df.m$TL, pch = 17, cex = 0.8, xlab = "", ylab = "", axes = FALSE, col = cmocean('haline')(4)[3], xlim = c(70, 175)) plot(mod.s, add=TRUE, link=FALSE, lwd=2, col = cmocean('haline')(4)[3]) axis(2, line = 8.2, las = 1) mtext(text = expression(paste(delta^34, "S")), side = 2, line = 10.5) abline(v = 85.06191, lty = "dashed", lwd = 2) legend(x = 70, y = 12, pch = c(15, 16, 17), col = c(cmocean('haline')(4)[1], cmocean('haline')(4)[2], cmocean('haline')(4)[3]), legend = c(expression(paste(delta^15, "N")), expression(paste(delta^13, "C")), expression(paste(delta^34, "S"))), cex = 0.8) # Liver par(mar = c(5, 12, 3, 1)) plot(df.l$d15NAIR ~ df.l$TL, subset = df.l$TL > 80, pch = 15, cex = 0.8, xlab = "Total length (cm)", ylab = "", axes = FALSE, col = cmocean('haline')(4)[1], xlim = c(70, 175)) axis(1, at = seq(70, 170, 10)) axis(2, las = 1) axis(3, at = c(70, t.save), labels = c("0","1", "", "2", "", "3", "", "4", "", "5", "", "6", "")) box() plot(mod.n2, add=TRUE, link=FALSE, lwd=2, col = cmocean('haline')(4)[1]) mtext(text = expression(paste(delta^15, "N")), side = 2, line = 2.5) par(new = TRUE) plot(df.l$d13CVPDB ~ df.l$TL, subset = df.l$TL > 80, pch = 16, cex = 0.8, xlab = "", ylab = "", axes = FALSE, col = cmocean('haline')(4)[2], xlim = c(70, 175)) plot(mod.c2, add=TRUE, link=FALSE, lwd=2, col = cmocean('haline')(4)[2]) axis(2, line = 4, las = 1) mtext(text = expression(paste(delta^13, "C")), side = 2, line = 6.5) par(new = TRUE) plot(df.l$d34SVCDT ~ df.l$TL, subset = df.l$TL > 80, pch = 17, cex = 0.8, xlab = "", ylab = "", axes = FALSE, col = cmocean('haline')(4)[3], xlim = c(70, 175)) plot(mod.s2, add=TRUE, link=FALSE, lwd=2, col = cmocean('haline')(4)[3]) axis(2, line = 8.2, las = 1) mtext(text = expression(paste(delta^34, "S")), side = 2, line = 10.5) abline(v = 85.06191, lty = "dashed", lwd = 2) legend(x = 70, y = 16, pch = c(15, 16, 17), col = c(cmocean('haline')(4)[1], cmocean('haline')(4)[2], cmocean('haline')(4)[3]), legend = c(expression(paste(delta^15, "N")), expression(paste(delta^13, "C")), expression(paste(delta^34, "S"))), cex = 0.8) dev.off()