--- title: "Breakpoint analyses" author: "Matthew Clements" date: "`r Sys.Date()`" output: html_document --- # File description ```{r} #This rmarkdown file are data analyses on crown-of-thorns competent brachiolaria larvae survival in response to low salinity, including levels that exist in GBR runoff plumes. ``` # Load packages ```{r} library(segmented) library(broom) library(brglm2) library(MASS) library(visreg) library(ResourceSelection) library(readxl) library(dplyr) library(strucchange) library(car) library(emmeans) library(knitr) library(ggplot2) library(gridExtra) library(grid) library(cowplot) library(lme4) library(MuMIn) library(ggplot2) library(gridExtra) ``` # ==================== # Breakpoint analyses # ==================== # Load data ```{r} #Piecewise and breakpoint analyses Data<-read_excel('[insert filepath]/Clements & Byrne 2026.xlsx', sheet = "Breakpoint analyses & SPCs", range = "A2:G147", col_names = c("Pre-settlement assay salinity exposure duration (days)", "Salinity_Treatment", "Startin n", "Normal", "Survival", "Normal (%)", "Survival (%)")) ``` ## Filter data by 'duration of 'Pre-settlement assay salinity exposure duration (days)' ```{r} Brach_2day_salinityexposure <- Data %>% filter(`Pre-settlement assay salinity exposure duration (days)` == "2") Brach_4day_salinityexposure <- Data %>% filter(`Pre-settlement assay salinity exposure duration (days)` == "4") Brach_7day_salinityexposure <- Data %>% filter(`Pre-settlement assay salinity exposure duration (days)` == "7") ``` ## Normal morphology and survival mean summary ### 2 day larval salinity exposure ```{r} # ===================== # Supplementary Table 2 # ===================== perc_normal_Brach_2day_salinityexposure <- Brach_2day_salinityexposure %>% filter(Salinity_Treatment %in% c(17, 19, 22, 25, 30, 34)) %>% group_by(Salinity_Treatment) %>% summarise( Mean = mean(`Normal (%)`, na.rm = TRUE), Std_Error = sd(`Normal (%)`, na.rm = TRUE) / sqrt(n()), count = n()) perc_normal_Brach_2day_salinityexposure perc_survival_Brach_2day_salinityexposure <- Brach_2day_salinityexposure %>% filter(Salinity_Treatment %in% c(17, 19, 22, 25, 30, 34)) %>% group_by(Salinity_Treatment) %>% summarise( Mean = mean(`Survival (%)`, na.rm = TRUE), Std_Error = sd(`Survival (%)`, na.rm = TRUE) / sqrt(n()), count = n()) perc_survival_Brach_2day_salinityexposure ``` ### 4 day larval salinity exposure ```{r} # ===================== # Supplementary Table 2 # ===================== perc_normal_Brach_4day_salinityexposure <- Brach_4day_salinityexposure %>% filter(Salinity_Treatment %in% c(17, 19, 22, 25, 30, 34)) %>% group_by(Salinity_Treatment) %>% summarise( Mean = mean(`Normal (%)`, na.rm = TRUE), Std_Error = sd(`Normal (%)`, na.rm = TRUE) / sqrt(n()), count = n()) perc_normal_Brach_4day_salinityexposure perc_survival_Brach_4day_salinityexposure <- Brach_4day_salinityexposure %>% filter(Salinity_Treatment %in% c(17, 19, 22, 25, 30, 34)) %>% group_by(Salinity_Treatment) %>% summarise( Mean = mean(`Survival (%)`, na.rm = TRUE), Std_Error = sd(`Survival (%)`, na.rm = TRUE) / sqrt(n()), count = n()) perc_survival_Brach_4day_salinityexposure ``` ### 7 day larval salinity exposure ```{r} # ===================== # Supplementary Table 2 # ===================== perc_normal_Brach_7day_salinityexposure <- Brach_7day_salinityexposure %>% filter(Salinity_Treatment %in% c(17, 19, 22, 25, 30, 34)) %>% group_by(Salinity_Treatment) %>% summarise( Mean = mean(`Normal (%)`, na.rm = TRUE), Std_Error = sd(`Normal (%)`, na.rm = TRUE) / sqrt(n()), count = n()) perc_normal_Brach_7day_salinityexposure perc_survival_Brach_7day_salinityexposure <- Brach_7day_salinityexposure %>% filter(Salinity_Treatment %in% c(17, 19, 22, 25, 30, 34)) %>% group_by(Salinity_Treatment) %>% summarise( Mean = mean(`Survival (%)`, na.rm = TRUE), Std_Error = sd(`Survival (%)`, na.rm = TRUE) / sqrt(n()), count = n()) perc_survival_Brach_7day_salinityexposure ``` ## Piecewise regressions and breakpoint analyses ### 2 day larval salinity exposure #### Normal larval morphology ```{r} #================================= # Model, pseudo r^2, and LS50/Sopt #================================= #Breakpoint detection Brach_2day_salinityexposure$norm_Proportion <- Brach_2day_salinityexposure$Normal / Brach_2day_salinityexposure$`Startin n` model_2day_salinityexposure <- segmented( glm(norm_Proportion ~ Salinity_Treatment, data = Brach_2day_salinityexposure, family = binomial(link = "logit"), method = "brglmFit"), seg.Z = ~ Salinity_Treatment, psi = 25 ) breakpoint_value <- model_2day_salinityexposure$psi[2] print(breakpoint_value) summary(model_2day_salinityexposure) # pseudo-rsqr logLik_model = logLik(model_2day_salinityexposure) null_model = glm(norm_Proportion ~ 1, data = Brach_2day_salinityexposure, family = binomial(link = "logit")) logLik_null = logLik(null_model) R2_McFadden = 1 - (logLik_model / logLik_null) print(R2_McFadden) # LS50 & Sopt90 new_data <- data.frame(Salinity_Treatment = seq(min(Data$Salinity_Treatment), max(Data$Salinity_Treatment), length.out = 1000)) new_data$predicted_prob <- predict(model_2day_salinityexposure, newdata = new_data, type = "response") # Predict probabilities using the segmented model target_probability <- 0.9 closest_index <- which.min(abs(new_data$predicted_prob - target_probability)) # Find the salinity treatment dose where the predicted probability is closest to the target probability dose_at_target_prob <- new_data$Salinity_Treatment[closest_index] print(dose_at_target_prob) target_probability <- 0.5 closest_index <- which.min(abs(new_data$predicted_prob - target_probability)) dose_at_target_prob <- new_data$Salinity_Treatment[closest_index] print(dose_at_target_prob) # =========== # Assumptions # =========== # 1. Linearity of the Logit visreg(model_2day_salinityexposure, "Salinity_Treatment", scale = "response", gg = TRUE) # 2. Goodness-of-Fit - Use Hosmer-Lemeshow test for goodness-of-fit hoslem.test(Brach_2day_salinityexposure$norm_Proportion, fitted(model_2day_salinityexposure), g = 10) # 3. Check for Overdispersion overdispersion <- sum(residuals(model_2day_salinityexposure, type = "pearson")^2) / df.residual(model_2day_salinityexposure) print(overdispersion) # 4. Diagnostic plots for residuals par(mfrow = c(2, 2)) plot(model_2day_salinityexposure) ``` ##### Create the prediction curve (breakpoint) and boxplot overlay ```{r} #======================= # Create the model curve #======================= # Generate a sequence of salinity treatments for predictions within the observed range salinity_range_Brach_2day_salinityexposure <- range(Brach_2day_salinityexposure$Salinity_Treatment, na.rm = TRUE) salinity_seq_Brach_2day_salinityexposure <- seq(from = salinity_range_Brach_2day_salinityexposure[1], to = salinity_range_Brach_2day_salinityexposure[2], by = 0.1) # Create a new data frame for predictions pred_data_Brach_2day_salinityexposure <- data.frame(Salinity_Treatment = salinity_seq_Brach_2day_salinityexposure) # Predict survival proportions using the model for 2 day salinity exposure data pred_data_Brach_2day_salinityexposure$Normal_Proportion_Predicted <- predict(model_2day_salinityexposure, newdata = pred_data_Brach_2day_salinityexposure, type = "response") # Create the prediction curve (breakpoint) breakpoint_curve_Brach_2day_salinityexposure <- ggplot(data = pred_data_Brach_2day_salinityexposure) + geom_line(aes(x = Salinity_Treatment, y = Normal_Proportion_Predicted), color = "red", size = 1) + geom_vline(xintercept = 19.06, linetype = "dashed", color = "grey", size = 1) + scale_x_continuous(breaks = c(17, 19, 22, 25, 30, 34)) + # Custom breaks on numeric x-axis scale_y_continuous(breaks = seq(0, 1, by = 0.1)) + # More tick marks on y-axis xlab("Salinity ‰") + ylab("Prob(Normal)") + ggtitle("2 d") + theme( plot.title = element_text(hjust = 0.5, size = 16, family = "Times New Roman", face = "bold"), panel.grid.major = element_line(color = "grey", size = 0.5), # Add major grid lines panel.grid.minor = element_line(color = "grey", size = 0.25), # Add minor grid lines panel.background = element_rect(fill = "white"), # Set background color to white axis.line = element_line(color = "black"), axis.text = element_text(size = 12, family = "Times New Roman"), axis.title.x = element_text(size = 16, family = "Times New Roman"), axis.title.y = element_text(size = 16, family = "Times New Roman") ) breakpoint_curve_Brach_2day_salinityexposure #========= # Box plot #========= # Assuming 'Salinity_Treatment' is the correct column that needs to be factored and converted to numeric Brach_2day_salinityexposure$Salinity_Treatment_Factor <- factor(Brach_2day_salinityexposure$Salinity_Treatment) Brach_2day_salinityexposure$Salinity_Treatment_Numeric <- as.numeric(as.character(Brach_2day_salinityexposure$Salinity_Treatment)) # Create a dataframe that will store the unique breaks and labels for the x-axis breaks_and_labels <- data.frame( Breaks = sort(unique(Brach_2day_salinityexposure$Salinity_Treatment_Numeric)), Labels = levels(Brach_2day_salinityexposure$Salinity_Treatment_Factor) ) # If there are more breaks than labels, this will truncate the breaks to match the number of labels breaks_and_labels <- breaks_and_labels[1:length(breaks_and_labels$Labels), ] # Plotting Boxplot_Brach_2day_salinityexposure <- ggplot(Brach_2day_salinityexposure, aes(x = Salinity_Treatment_Numeric, y = `Normal (%)`, group = Salinity_Treatment_Factor)) + geom_boxplot(aes(fill = Salinity_Treatment_Factor), show.legend = FALSE, outlier.shape = TRUE, width = 1.5) + scale_x_continuous(name = "Salinity (‰)", breaks = breaks_and_labels$Breaks, labels = breaks_and_labels$Labels) + scale_y_continuous(name = "Normal (%)", limits = c(0, 100)) + scale_fill_manual(values = rep("#bad6f9", length(unique(Brach_2day_salinityexposure$Salinity_Treatment_Factor)))) + ggtitle("2 d") + theme(plot.title = element_text(hjust = 0.5, size = 16, family = "Times New Roman", face = "bold"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), axis.text = element_text(size = 12, family = "Times New Roman"), axis.title = element_text(size = 16, family = "Times New Roman")) # Display the plot print(Boxplot_Brach_2day_salinityexposure) #==================================== # Overlay the model onto the box plot #==================================== Brach_breakpoint_curves_Brach_2day_salinityexposure <- plot_grid(Brach_2day_salinityexposure, breakpoint_curve_Brach_2day_salinityexposure, ncol = 2, nrow = 1) Brach_bps <- grid.arrange(Brach_breakpoint_curves_Brach_2day_salinityexposure) # Ensure that pred_data Salinity_Treatment is numeric and matches the Datas numeric scale pred_data_Brach_2day_salinityexposure$Salinity_Treatment_Numeric <- as.numeric(as.character(pred_data_Brach_2day_salinityexposure$Salinity_Treatment)) Figure_3a <- ggplot(Brach_2day_salinityexposure, aes(x = Salinity_Treatment_Numeric, y = `Normal (%)`)) + geom_boxplot(aes(fill = factor(Salinity_Treatment_Numeric)), show.legend = FALSE, outlier.shape = TRUE, width = 1.5) + scale_fill_manual(values = rep("#f0f0f0", length(unique(Brach_2day_salinityexposure$Salinity_Treatment_Numeric)))) + scale_x_continuous(name = "", breaks = sort(unique(Brach_2day_salinityexposure$Salinity_Treatment_Numeric)), labels = sort(unique(Brach_2day_salinityexposure$Salinity_Treatment_Numeric))) + scale_y_continuous( name = "Normal (%)", limits = c(0, 100), sec.axis = sec_axis(~ . / 100, name = "") # Define secondary axis here ) + geom_line(data = pred_data_Brach_2day_salinityexposure, aes(x = Salinity_Treatment_Numeric, y = Normal_Proportion_Predicted * 100), color = "black", size = 0.5) + geom_vline(xintercept = 19.06, linetype = "dashed", color = "black", size = 0.7) + geom_segment(aes(x = 27.9, y = 5, xend = 27.9, yend = 0), arrow = arrow(type = "closed", length = unit(0.2, "inches")), colour = "#525252", size = 0.5) + # Arrow at 27.9 pointing downwards geom_segment(aes(x = 19, y = 5, xend = 19, yend = 0), arrow = arrow(type = "closed", length = unit(0.2, "inches")), colour = "#bdbdbd", size = 0.5) + # Arrow at 19 pointing downwards ggtitle("2 days in salinity", "(a)") + theme( plot.title = element_text(hjust = 0.5, size = 30, family = "Times New Roman", face = "bold"), plot.subtitle = element_text(hjust = 0, size = 30, family = "Times New Roman", face = "bold"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), axis.text = element_text(size = 30, family = "Times New Roman"), axis.title = element_text(size = 30, family = "Times New Roman"), axis.title.y.right = element_text(color = "black") ) # Print the plot with arrows pointing downwards print(Figure_3a) greyscale_colors <- c("#f0f0f0", "#d9d9d9", "#bdbdbd", "#969696", "#737373", "#525252", "#252525") ``` #### Survival ```{r} #================================= # Model, pseudo r^2, and LS50/Sopt #================================= #Breakpoint detection Brach_2day_salinityexposure$Surv_prop <- Brach_2day_salinityexposure$Survival / Brach_2day_salinityexposure$`Startin n` model_2day_salinityexposure <- segmented( glm(Surv_prop ~ Salinity_Treatment, data = Brach_2day_salinityexposure, family = binomial(link = "logit"), method = "brglmFit"), seg.Z = ~ Salinity_Treatment, psi = 25 ) breakpoint_value <- model_2day_salinityexposure$psi[2] print(breakpoint_value) summary(model_2day_salinityexposure) # pseudo-rsqr logLik_model = logLik(model_2day_salinityexposure) null_model = glm(Surv_prop ~ 1, data = Brach_2day_salinityexposure, family = binomial(link = "logit")) logLik_null = logLik(null_model) R2_McFadden = 1 - (logLik_model / logLik_null) print(R2_McFadden) # LS50 & Sopt90 new_data <- data.frame(Salinity_Treatment = seq(min(Data$Salinity_Treatment), max(Data$Salinity_Treatment), length.out = 1000)) new_data$predicted_prob <- predict(model_2day_salinityexposure, newdata = new_data, type = "response") # Predict probabilities using the segmented model target_probability <- 0.9 closest_index <- which.min(abs(new_data$predicted_prob - target_probability)) # Find the salinity treatment dose where the predicted probability is closest to the target probability dose_at_target_prob <- new_data$Salinity_Treatment[closest_index] print(dose_at_target_prob) target_probability <- 0.5 closest_index <- which.min(abs(new_data$predicted_prob - target_probability)) dose_at_target_prob <- new_data$Salinity_Treatment[closest_index] print(dose_at_target_prob) # =========== # Assumptions # =========== # 1. Linearity of the Logit visreg(model_2day_salinityexposure, "Salinity_Treatment", scale = "response", gg = TRUE) # 2. Goodness-of-Fit - Use Hosmer-Lemeshow test for goodness-of-fit hoslem.test(Brach_2day_salinityexposure$Surv_prop, fitted(model_2day_salinityexposure), g = 10) # 3. Check for Overdispersion overdispersion <- sum(residuals(model_2day_salinityexposure, type = "pearson")^2) / df.residual(model_2day_salinityexposure) print(overdispersion) # 4. Diagnostic plots for residuals par(mfrow = c(2, 2)) plot(model_2day_salinityexposure) ``` ##### Create the prediction curve (breakpoint) and boxplot overlay ```{r} #======================= # Create the model curve #======================= # Generate a sequence of salinity treatments for predictions within the observed range salinity_range_Brach_2day_salinityexposure <- range(Brach_2day_salinityexposure$Salinity_Treatment, na.rm = TRUE) salinity_seq_Brach_2day_salinityexposure <- seq(from = salinity_range_Brach_2day_salinityexposure[1], to = salinity_range_Brach_2day_salinityexposure[2], by = 0.1) # Create a new data frame for predictions pred_data_Brach_2day_salinityexposure <- data.frame(Salinity_Treatment = salinity_seq_Brach_2day_salinityexposure) # Predict survival proportions using the model for 2 day salinity exposure data pred_data_Brach_2day_salinityexposure$Surv_Proportion_Predicted <- predict(model_2day_salinityexposure, newdata = pred_data_Brach_2day_salinityexposure, type = "response") # Create the prediction curve (breakpoint) breakpoint_curve_Brach_2day_salinityexposure <- ggplot(data = pred_data_Brach_2day_salinityexposure) + geom_line(aes(x = Salinity_Treatment, y = Surv_Proportion_Predicted), color = "red", size = 1) + geom_vline(xintercept = 19.06, linetype = "dashed", color = "grey", size = 1) + scale_x_continuous(breaks = c(17, 19, 22, 25, 30, 34)) + # Custom breaks on numeric x-axis scale_y_continuous(breaks = seq(0, 1, by = 0.1)) + # More tick marks on y-axis xlab("Salinity ‰") + ylab("Prob(Normal)") + ggtitle("2 d") + theme( plot.title = element_text(hjust = 0.5, size = 16, family = "Times New Roman", face = "bold"), panel.grid.major = element_line(color = "grey", size = 0.5), # Add major grid lines panel.grid.minor = element_line(color = "grey", size = 0.25), # Add minor grid lines panel.background = element_rect(fill = "white"), # Set background color to white axis.line = element_line(color = "black"), axis.text = element_text(size = 12, family = "Times New Roman"), axis.title.x = element_text(size = 16, family = "Times New Roman"), axis.title.y = element_text(size = 16, family = "Times New Roman") ) breakpoint_curve_Brach_2day_salinityexposure #========= # Box plot #========= # Assuming 'Salinity_Treatment' is the correct column that needs to be factored and converted to numeric Brach_2day_salinityexposure$Salinity_Treatment_Factor <- factor(Brach_2day_salinityexposure$Salinity_Treatment) Brach_2day_salinityexposure$Salinity_Treatment_Numeric <- as.numeric(as.character(Brach_2day_salinityexposure$Salinity_Treatment)) # Create a dataframe that will store the unique breaks and labels for the x-axis breaks_and_labels <- data.frame( Breaks = sort(unique(Brach_2day_salinityexposure$Salinity_Treatment_Numeric)), Labels = levels(Brach_2day_salinityexposure$Salinity_Treatment_Factor) ) # If there are more breaks than labels, this will truncate the breaks to match the number of labels breaks_and_labels <- breaks_and_labels[1:length(breaks_and_labels$Labels), ] # Plotting Boxplot_Brach_2day_salinityexposure <- ggplot(Brach_2day_salinityexposure, aes(x = Salinity_Treatment_Numeric, y = `Survival (%)`, group = Salinity_Treatment_Factor)) + geom_boxplot(aes(fill = Salinity_Treatment_Factor), show.legend = FALSE, outlier.shape = TRUE, width = 1.5) + scale_x_continuous(name = "Salinity (‰)", breaks = breaks_and_labels$Breaks, labels = breaks_and_labels$Labels) + scale_y_continuous(name = "Survival (%)", limits = c(0, 100)) + scale_fill_manual(values = rep("#bad6f9", length(unique(Brach_2day_salinityexposure$Salinity_Treatment_Factor)))) + ggtitle("2 d") + theme(plot.title = element_text(hjust = 0.5, size = 16, family = "Times New Roman", face = "bold"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), axis.text = element_text(size = 12, family = "Times New Roman"), axis.title = element_text(size = 16, family = "Times New Roman")) # Display the plot print(Boxplot_Brach_2day_salinityexposure) #==================================== # Overlay the model onto the box plot #==================================== Brach_breakpoint_curves_Brach_2day_salinityexposure <- plot_grid(Brach_2day_salinityexposure, breakpoint_curve_Brach_2day_salinityexposure, ncol = 2, nrow = 1) Brach_bps <- grid.arrange(Brach_breakpoint_curves_Brach_2day_salinityexposure) # Ensure that pred_data Salinity_Treatment is numeric and matches the Datas numeric scale pred_data_Brach_2day_salinityexposure$Salinity_Treatment_Numeric <- as.numeric(as.character(pred_data_Brach_2day_salinityexposure$Salinity_Treatment)) Figure_3d <- ggplot(Brach_2day_salinityexposure, aes(x = Salinity_Treatment_Numeric, y = `Survival (%)`)) + geom_boxplot(aes(fill = factor(Salinity_Treatment_Numeric)), show.legend = FALSE, outlier.shape = TRUE, width = 1.5) + scale_fill_manual(values = rep("#f0f0f0", length(unique(Brach_2day_salinityexposure$Salinity_Treatment_Numeric)))) + scale_x_continuous(name = "", breaks = sort(unique(Brach_2day_salinityexposure$Salinity_Treatment_Numeric)), labels = sort(unique(Brach_2day_salinityexposure$Salinity_Treatment_Numeric))) + scale_y_continuous( name = "Survival (%)", limits = c(0, 100), sec.axis = sec_axis(~ . / 100, name = "") ) + geom_line(data = pred_data_Brach_2day_salinityexposure, aes(x = Salinity_Treatment_Numeric, y = Surv_Proportion_Predicted * 100), color = "black", size = 0.5) + geom_vline(xintercept = 19.2, linetype = "dashed", color = "black", size = 0.7) + geom_segment(aes(x = 22, y = 5, xend = 22, yend = 0), arrow = arrow(type = "closed", length = unit(0.2, "inches")), colour = "#525252", size = 0.5) + geom_segment(aes(x = 18, y = 5, xend = 18, yend = 0), arrow = arrow(type = "closed", length = unit(0.2, "inches")), colour = "#bdbdbd", size = 0.5) + ggtitle("", "(d)") + theme( plot.title = element_text(hjust = 0.5, size = 30, family = "Times New Roman", face = "bold"), plot.subtitle = element_text(hjust = 0, size = 30, family = "Times New Roman", face = "bold"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), axis.text = element_text(size = 30, family = "Times New Roman"), axis.title = element_text(size = 30, family = "Times New Roman"), axis.title.y.right = element_text(color = "black") ) # Print the plot with arrows pointing downwards print(Figure_3d) greyscale_colors <- c("#f0f0f0", "#d9d9d9", "#bdbdbd", "#969696", "#737373", "#525252", "#252525") ``` ### 4 day larval salinity exposure #### Normal larval morphology ```{r} #================================= # Model, pseudo r^2, and LS50/Sopt #================================= #Breakpoint detection Brach_4day_salinityexposure$norm_Proportion <- Brach_4day_salinityexposure$Normal / Brach_4day_salinityexposure$`Startin n` model_4day_salinityexposure <- segmented( glm(norm_Proportion ~ Salinity_Treatment, data = Brach_4day_salinityexposure, family = binomial(link = "logit"), method = "brglmFit"), seg.Z = ~ Salinity_Treatment, psi = 25 ) breakpoint_value <- model_4day_salinityexposure$psi[2] print(breakpoint_value) summary(model_4day_salinityexposure) # pseudo-rsqr logLik_model = logLik(model_4day_salinityexposure) null_model = glm(norm_Proportion ~ 1, data = Brach_4day_salinityexposure, family = binomial(link = "logit")) logLik_null = logLik(null_model) R2_McFadden = 1 - (logLik_model / logLik_null) print(R2_McFadden) # LS50 & Sopt90 new_data <- data.frame(Salinity_Treatment = seq(min(Data$Salinity_Treatment), max(Data$Salinity_Treatment), length.out = 1000)) new_data$predicted_prob <- predict(model_4day_salinityexposure, newdata = new_data, type = "response") # Predict probabilities using the segmented model target_probability <- 0.9 closest_index <- which.min(abs(new_data$predicted_prob - target_probability)) # Find the salinity treatment dose where the predicted probability is closest to the target probability dose_at_target_prob <- new_data$Salinity_Treatment[closest_index] print(dose_at_target_prob) target_probability <- 0.5 closest_index <- which.min(abs(new_data$predicted_prob - target_probability)) dose_at_target_prob <- new_data$Salinity_Treatment[closest_index] print(dose_at_target_prob) # =========== # Assumptions # =========== # 1. Linearity of the Logit visreg(model_4day_salinityexposure, "Salinity_Treatment", scale = "response", gg = TRUE) # 2. Goodness-of-Fit - Use Hosmer-Lemeshow test for goodness-of-fit hoslem.test(Brach_4day_salinityexposure$norm_Proportion, fitted(model_4day_salinityexposure), g = 10) # 3. Check for Overdispersion overdispersion <- sum(residuals(model_4day_salinityexposure, type = "pearson")^2) / df.residual(model_4day_salinityexposure) print(overdispersion) # 4. Diagnostic plots for residuals par(mfrow = c(2, 2)) plot(model_4day_salinityexposure) ``` ##### Create the prediction curve (breakpoint) and boxplot overlay ```{r} #======================= # Create the model curve #======================= # Generate a sequence of salinity treatments for predictions within the observed range salinity_range_Brach_4day_salinityexposure <- range(Brach_4day_salinityexposure$Salinity_Treatment, na.rm = TRUE) salinity_seq_Brach_4day_salinityexposure <- seq(from = salinity_range_Brach_4day_salinityexposure[1], to = salinity_range_Brach_4day_salinityexposure[2], by = 0.1) # Create a new data frame for predictions pred_data_Brach_4day_salinityexposure <- data.frame(Salinity_Treatment = salinity_seq_Brach_4day_salinityexposure) # Predict survival proportions using the model for 4 day salinity exposure data pred_data_Brach_4day_salinityexposure$Normal_Proportion_Predicted <- predict(model_4day_salinityexposure, newdata = pred_data_Brach_4day_salinityexposure, type = "response") # Create the prediction curve (breakpoint) breakpoint_curve_Brach_4day_salinityexposure <- ggplot(data = pred_data_Brach_4day_salinityexposure) + geom_line(aes(x = Salinity_Treatment, y = Normal_Proportion_Predicted), color = "red", size = 1) + geom_vline(xintercept = 22.6, linetype = "dashed", color = "grey", size = 1) + scale_x_continuous(breaks = c(17, 19, 22, 25, 30, 34)) + # Custom breaks on numeric x-axis scale_y_continuous(breaks = seq(0, 1, by = 0.1)) + # More tick marks on y-axis xlab("Salinity ‰") + ylab("Prob(Normal)") + ggtitle("4 d") + theme( plot.title = element_text(hjust = 0.5, size = 16, family = "Times New Roman", face = "bold"), panel.grid.major = element_line(color = "grey", size = 0.5), # Add major grid lines panel.grid.minor = element_line(color = "grey", size = 0.25), # Add minor grid lines panel.background = element_rect(fill = "white"), # Set background color to white axis.line = element_line(color = "black"), axis.text = element_text(size = 12, family = "Times New Roman"), axis.title.x = element_text(size = 16, family = "Times New Roman"), axis.title.y = element_text(size = 16, family = "Times New Roman") ) breakpoint_curve_Brach_4day_salinityexposure #========= # Box plot #========= # Assuming 'Salinity_Treatment' is the correct column that needs to be factored and converted to numeric Brach_4day_salinityexposure$Salinity_Treatment_Factor <- factor(Brach_4day_salinityexposure$Salinity_Treatment) Brach_4day_salinityexposure$Salinity_Treatment_Numeric <- as.numeric(as.character(Brach_4day_salinityexposure$Salinity_Treatment)) # Create a dataframe that will store the unique breaks and labels for the x-axis breaks_and_labels <- data.frame( Breaks = sort(unique(Brach_4day_salinityexposure$Salinity_Treatment_Numeric)), Labels = levels(Brach_4day_salinityexposure$Salinity_Treatment_Factor) ) # If there are more breaks than labels, this will truncate the breaks to match the number of labels breaks_and_labels <- breaks_and_labels[1:length(breaks_and_labels$Labels), ] # Plotting Boxplot_Brach_4day_salinityexposure <- ggplot(Brach_4day_salinityexposure, aes(x = Salinity_Treatment_Numeric, y = `Normal (%)`, group = Salinity_Treatment_Factor)) + geom_boxplot(aes(fill = Salinity_Treatment_Factor), show.legend = FALSE, outlier.shape = TRUE, width = 1.5) + scale_x_continuous(name = "Salinity (‰)", breaks = breaks_and_labels$Breaks, labels = breaks_and_labels$Labels) + scale_y_continuous(name = "Normal (%)", limits = c(0, 100)) + scale_fill_manual(values = rep("#bad6f9", length(unique(Brach_4day_salinityexposure$Salinity_Treatment_Factor)))) + ggtitle("4 d") + theme(plot.title = element_text(hjust = 0.5, size = 16, family = "Times New Roman", face = "bold"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), axis.text = element_text(size = 12, family = "Times New Roman"), axis.title = element_text(size = 16, family = "Times New Roman")) # Display the plot print(Boxplot_Brach_4day_salinityexposure) #==================================== # Overlay the model onto the box plot #==================================== Brach_breakpoint_curves_Brach_4day_salinityexposure <- plot_grid(Brach_4day_salinityexposure, breakpoint_curve_Brach_4day_salinityexposure, ncol = 2, nrow = 1) Brach_bps <- grid.arrange(Brach_breakpoint_curves_Brach_4day_salinityexposure) # Ensure that pred_data Salinity_Treatment is numeric and matches the Datas numeric scale pred_data_Brach_4day_salinityexposure$Salinity_Treatment_Numeric <- as.numeric(as.character(pred_data_Brach_4day_salinityexposure$Salinity_Treatment)) Figure_3b <- ggplot(Brach_4day_salinityexposure, aes(x = Salinity_Treatment_Numeric, y = `Normal (%)`)) + geom_boxplot(aes(fill = factor(Salinity_Treatment_Numeric)), show.legend = FALSE, outlier.shape = TRUE, width = 1.5) + scale_fill_manual(values = rep("#f0f0f0", length(unique(Brach_4day_salinityexposure$Salinity_Treatment_Numeric)))) + scale_x_continuous(name = "", breaks = sort(unique(Brach_4day_salinityexposure$Salinity_Treatment_Numeric)), labels = sort(unique(Brach_4day_salinityexposure$Salinity_Treatment_Numeric))) + scale_y_continuous( name = "", limits = c(0, 100), sec.axis = sec_axis(~ . / 100, name = "") # Define secondary axis here ) + geom_line(data = pred_data_Brach_4day_salinityexposure, aes(x = Salinity_Treatment_Numeric, y = Normal_Proportion_Predicted * 100), color = "black", size = 0.5) + geom_vline(xintercept = 22.6, linetype = "dashed", color = "black", size = 0.7) + geom_segment(aes(x = 28.6, y = 5, xend = 28.6, yend = 0), arrow = arrow(type = "closed", length = unit(0.2, "inches")), colour = "#525252", size = 0.5) + # Arrow at 27.9 pointing downwards geom_segment(aes(x = 22.4, y = 5, xend = 22.4, yend = 0), arrow = arrow(type = "closed", length = unit(0.2, "inches")), colour = "#bdbdbd", size = 0.5) + # Arrow at 19 pointing downwards ggtitle("4 days in salinity", "(b)") + theme( plot.title = element_text(hjust = 0.5, size = 30, family = "Times New Roman", face = "bold"), plot.subtitle = element_text(hjust = 0, size = 30, family = "Times New Roman", face = "bold"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), axis.text = element_text(size = 30, family = "Times New Roman"), axis.title = element_text(size = 30, family = "Times New Roman"), axis.title.y.right = element_text(color = "black") ) # Print the plot with arrows pointing downwards print(Figure_3b) greyscale_colors <- c("#f0f0f0", "#d9d9d9", "#bdbdbd", "#969696", "#737373", "#525252", "#252525") ``` #### Survival ```{r} #================================= # Model, pseudo r^2, and LS50/Sopt #================================= #Breakpoint detection Brach_4day_salinityexposure$Surv_prop <- Brach_4day_salinityexposure$Survival / Brach_4day_salinityexposure$`Startin n` model_4day_salinityexposure <- segmented( glm(Surv_prop ~ Salinity_Treatment, data = Brach_4day_salinityexposure, family = binomial(link = "logit")), seg.Z = ~ Salinity_Treatment, psi = 25 ) breakpoint_value <- model_4day_salinityexposure$psi[2] print(breakpoint_value) summary(model_4day_salinityexposure) # pseudo-rsqr logLik_model = logLik(model_4day_salinityexposure) null_model = glm(Surv_prop ~ 1, data = Brach_4day_salinityexposure, family = binomial(link = "logit")) logLik_null = logLik(null_model) R2_McFadden = 1 - (logLik_model / logLik_null) print(R2_McFadden) # LS50 & Sopt90 new_data <- data.frame(Salinity_Treatment = seq(min(Data$Salinity_Treatment), max(Data$Salinity_Treatment), length.out = 1000)) new_data$predicted_prob <- predict(model_4day_salinityexposure, newdata = new_data, type = "response") # Predict probabilities using the segmented model target_probability <- 0.9 closest_index <- which.min(abs(new_data$predicted_prob - target_probability)) # Find the salinity treatment dose where the predicted probability is closest to the target probability dose_at_target_prob <- new_data$Salinity_Treatment[closest_index] print(dose_at_target_prob) target_probability <- 0.5 closest_index <- which.min(abs(new_data$predicted_prob - target_probability)) dose_at_target_prob <- new_data$Salinity_Treatment[closest_index] print(dose_at_target_prob) # =========== # Assumptions # =========== # 1. Linearity of the Logit visreg(model_4day_salinityexposure, "Salinity_Treatment", scale = "response", gg = TRUE) # 2. Goodness-of-Fit - Use Hosmer-Lemeshow test for goodness-of-fit hoslem.test(Brach_4day_salinityexposure$Surv_prop, fitted(model_4day_salinityexposure), g = 10) # 3. Check for Overdispersion overdispersion <- sum(residuals(model_4day_salinityexposure, type = "pearson")^2) / df.residual(model_4day_salinityexposure) print(overdispersion) # 4. Diagnostic plots for residuals par(mfrow = c(2, 2)) plot(model_4day_salinityexposure) ``` ##### Create the prediction curve (breakpoint) and boxplot overlay ```{r} #======================= # Create the model curve #======================= # Generate a sequence of salinity treatments for predictions within the observed range salinity_range_Brach_4day_salinityexposure <- range(Brach_4day_salinityexposure$Salinity_Treatment, na.rm = TRUE) salinity_seq_Brach_4day_salinityexposure <- seq(from = salinity_range_Brach_4day_salinityexposure[1], to = salinity_range_Brach_4day_salinityexposure[2], by = 0.1) # Create a new data frame for predictions pred_data_Brach_4day_salinityexposure <- data.frame(Salinity_Treatment = salinity_seq_Brach_4day_salinityexposure) # Predict survival proportions using the model for 4 day salinity exposure data pred_data_Brach_4day_salinityexposure$Surv_Proportion_Predicted <- predict(model_4day_salinityexposure, newdata = pred_data_Brach_4day_salinityexposure, type = "response") # Create the prediction curve (breakpoint) breakpoint_curve_Brach_4day_salinityexposure <- ggplot(data = pred_data_Brach_4day_salinityexposure) + geom_line(aes(x = Salinity_Treatment, y = Surv_Proportion_Predicted), color = "red", size = 1) + geom_vline(xintercept = 19.8, linetype = "dashed", color = "grey", size = 1) + scale_x_continuous(breaks = c(17, 19, 22, 25, 30, 34)) + # Custom breaks on numeric x-axis scale_y_continuous(breaks = seq(0, 1, by = 0.1)) + # More tick marks on y-axis xlab("Salinity ‰") + ylab("Prob(Surv)") + ggtitle("4 d") + theme( plot.title = element_text(hjust = 0.5, size = 16, family = "Times New Roman", face = "bold"), panel.grid.major = element_line(color = "grey", size = 0.5), # Add major grid lines panel.grid.minor = element_line(color = "grey", size = 0.25), # Add minor grid lines panel.background = element_rect(fill = "white"), # Set background color to white axis.line = element_line(color = "black"), axis.text = element_text(size = 12, family = "Times New Roman"), axis.title.x = element_text(size = 16, family = "Times New Roman"), axis.title.y = element_text(size = 16, family = "Times New Roman") ) breakpoint_curve_Brach_4day_salinityexposure #========= # Box plot #========= # Assuming 'Salinity_Treatment' is the correct column that needs to be factored and converted to numeric Brach_4day_salinityexposure$Salinity_Treatment_Factor <- factor(Brach_4day_salinityexposure$Salinity_Treatment) Brach_4day_salinityexposure$Salinity_Treatment_Numeric <- as.numeric(as.character(Brach_4day_salinityexposure$Salinity_Treatment)) # Create a dataframe that will store the unique breaks and labels for the x-axis breaks_and_labels <- data.frame( Breaks = sort(unique(Brach_4day_salinityexposure$Salinity_Treatment_Numeric)), Labels = levels(Brach_4day_salinityexposure$Salinity_Treatment_Factor) ) # If there are more breaks than labels, this will truncate the breaks to match the number of labels breaks_and_labels <- breaks_and_labels[1:length(breaks_and_labels$Labels), ] # Plotting Boxplot_Brach_4day_salinityexposure <- ggplot(Brach_4day_salinityexposure, aes(x = Salinity_Treatment_Numeric, y = `Survival (%)`, group = Salinity_Treatment_Factor)) + geom_boxplot(aes(fill = Salinity_Treatment_Factor), show.legend = FALSE, outlier.shape = TRUE, width = 1.5) + scale_x_continuous(name = "Salinity (‰)", breaks = breaks_and_labels$Breaks, labels = breaks_and_labels$Labels) + scale_y_continuous(name = "Survival (%)", limits = c(0, 100)) + scale_fill_manual(values = rep("#bad6f9", length(unique(Brach_4day_salinityexposure$Salinity_Treatment_Factor)))) + ggtitle("4 d") + theme(plot.title = element_text(hjust = 0.5, size = 16, family = "Times New Roman", face = "bold"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), axis.text = element_text(size = 12, family = "Times New Roman"), axis.title = element_text(size = 16, family = "Times New Roman")) # Display the plot print(Boxplot_Brach_4day_salinityexposure) #==================================== # Overlay the model onto the box plot #==================================== Brach_breakpoint_curves_Brach_4day_salinityexposure <- plot_grid(Brach_4day_salinityexposure, breakpoint_curve_Brach_4day_salinityexposure, ncol = 2, nrow = 1) Brach_bps <- grid.arrange(Brach_breakpoint_curves_Brach_4day_salinityexposure) # Ensure that pred_data Salinity_Treatment is numeric and matches the Datas numeric scale pred_data_Brach_4day_salinityexposure$Salinity_Treatment_Numeric <- as.numeric(as.character(pred_data_Brach_4day_salinityexposure$Salinity_Treatment)) Figure_3e <- ggplot(Brach_4day_salinityexposure, aes(x = Salinity_Treatment_Numeric, y = `Survival (%)`)) + geom_boxplot(aes(fill = factor(Salinity_Treatment_Numeric)), show.legend = FALSE, outlier.shape = TRUE, width = 1.5) + scale_fill_manual(values = rep("#f0f0f0", length(unique(Brach_4day_salinityexposure$Salinity_Treatment_Numeric)))) + scale_x_continuous(name = "", breaks = sort(unique(Brach_4day_salinityexposure$Salinity_Treatment_Numeric)), labels = sort(unique(Brach_4day_salinityexposure$Salinity_Treatment_Numeric))) + scale_y_continuous( name = "", limits = c(0, 100), sec.axis = sec_axis(~ . / 100, name = "") ) + geom_line(data = pred_data_Brach_4day_salinityexposure, aes(x = Salinity_Treatment_Numeric, y = Surv_Proportion_Predicted * 100), color = "black", size = 0.5) + geom_vline(xintercept = 19.8, linetype = "dashed", color = "black", size = 0.7) + geom_segment(aes(x = 23.4, y = 5, xend = 23.4, yend = 0), arrow = arrow(type = "closed", length = unit(0.2, "inches")), colour = "#525252", size = 0.5) + geom_segment(aes(x = 17.8, y = 5, xend = 17.8, yend = 0), arrow = arrow(type = "closed", length = unit(0.2, "inches")), colour = "#bdbdbd", size = 0.5) + ggtitle("", "(e)") + theme( plot.title = element_text(hjust = 0.5, size = 30, family = "Times New Roman", face = "bold"), plot.subtitle = element_text(hjust = 0, size = 30, family = "Times New Roman", face = "bold"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), axis.text = element_text(size = 30, family = "Times New Roman"), axis.title = element_text(size = 30, family = "Times New Roman"), axis.title.y.right = element_text(color = "black") ) # Print the plot with arrows pointing downwards print(Figure_3e) greyscale_colors <- c("#f0f0f0", "#d9d9d9", "#bdbdbd", "#969696", "#737373", "#525252", "#252525") ``` ### 7 day larval salinity exposure #### Normal larval morphology ```{r} #================================= # Model, pseudo r^2, and LS50/Sopt #================================= #Breakpoint detection Brach_7day_salinityexposure$norm_Proportion <- Brach_7day_salinityexposure$Normal / Brach_7day_salinityexposure$`Startin n` model_7day_salinityexposure <- segmented( glm(norm_Proportion ~ Salinity_Treatment, data = Brach_7day_salinityexposure, family = binomial(link = "logit"), method = "brglmFit"), seg.Z = ~ Salinity_Treatment, psi = 25 ) breakpoint_value <- model_7day_salinityexposure$psi[2] print(breakpoint_value) summary(model_7day_salinityexposure) # pseudo-rsqr logLik_model = logLik(model_7day_salinityexposure) null_model = glm(norm_Proportion ~ 1, data = Brach_7day_salinityexposure, family = binomial(link = "logit")) logLik_null = logLik(null_model) R2_McFadden = 1 - (logLik_model / logLik_null) print(R2_McFadden) # LS50 & Sopt90 new_data <- data.frame(Salinity_Treatment = seq(min(Data$Salinity_Treatment), max(Data$Salinity_Treatment), length.out = 1000)) new_data$predicted_prob <- predict(model_7day_salinityexposure, newdata = new_data, type = "response") # Predict probabilities using the segmented model target_probability <- 0.9 closest_index <- which.min(abs(new_data$predicted_prob - target_probability)) # Find the salinity treatment dose where the predicted probability is closest to the target probability dose_at_target_prob <- new_data$Salinity_Treatment[closest_index] print(dose_at_target_prob) target_probability <- 0.5 closest_index <- which.min(abs(new_data$predicted_prob - target_probability)) dose_at_target_prob <- new_data$Salinity_Treatment[closest_index] print(dose_at_target_prob) # =========== # Assumptions # =========== # 1. Linearity of the Logit visreg(model_7day_salinityexposure, "Salinity_Treatment", scale = "response", gg = TRUE) # 2. Goodness-of-Fit - Use Hosmer-Lemeshow test for goodness-of-fit hoslem.test(Brach_7day_salinityexposure$norm_Proportion, fitted(model_7day_salinityexposure), g = 10) # 3. Check for Overdispersion overdispersion <- sum(residuals(model_7day_salinityexposure, type = "pearson")^2) / df.residual(model_7day_salinityexposure) print(overdispersion) # 4. Diagnostic plots for residuals par(mfrow = c(2, 2)) plot(model_7day_salinityexposure) ``` ##### Create the prediction curve (breakpoint) and boxplot overlay ```{r} #======================= # Create the model curve #======================= # Generate a sequence of salinity treatments for predictions within the observed range salinity_range_Brach_7day_salinityexposure <- range(Brach_7day_salinityexposure$Salinity_Treatment, na.rm = TRUE) salinity_seq_Brach_7day_salinityexposure <- seq(from = salinity_range_Brach_7day_salinityexposure[1], to = salinity_range_Brach_7day_salinityexposure[2], by = 0.1) # Create a new data frame for predictions pred_data_Brach_7day_salinityexposure <- data.frame(Salinity_Treatment = salinity_seq_Brach_7day_salinityexposure) # Predict survival proportions using the model for 4 day salinity exposure data pred_data_Brach_7day_salinityexposure$Normal_Proportion_Predicted <- predict(model_7day_salinityexposure, newdata = pred_data_Brach_7day_salinityexposure, type = "response") # Create the prediction curve (breakpoint) breakpoint_curve_Brach_7day_salinityexposure <- ggplot(data = pred_data_Brach_7day_salinityexposure) + geom_line(aes(x = Salinity_Treatment, y = Normal_Proportion_Predicted), color = "red", size = 1) + geom_vline(xintercept = 22.2, linetype = "dashed", color = "grey", size = 1) + scale_x_continuous(breaks = c(17, 19, 22, 25, 30, 34)) + # Custom breaks on numeric x-axis scale_y_continuous(breaks = seq(0, 1, by = 0.1)) + # More tick marks on y-axis xlab("Salinity ‰") + ylab("Prob(Normal)") + ggtitle("7 d") + theme( plot.title = element_text(hjust = 0.5, size = 16, family = "Times New Roman", face = "bold"), panel.grid.major = element_line(color = "grey", size = 0.5), # Add major grid lines panel.grid.minor = element_line(color = "grey", size = 0.25), # Add minor grid lines panel.background = element_rect(fill = "white"), # Set background color to white axis.line = element_line(color = "black"), axis.text = element_text(size = 12, family = "Times New Roman"), axis.title.x = element_text(size = 16, family = "Times New Roman"), axis.title.y = element_text(size = 16, family = "Times New Roman") ) breakpoint_curve_Brach_7day_salinityexposure #========= # Box plot #========= # Assuming 'Salinity_Treatment' is the correct column that needs to be factored and converted to numeric Brach_7day_salinityexposure$Salinity_Treatment_Factor <- factor(Brach_7day_salinityexposure$Salinity_Treatment) Brach_7day_salinityexposure$Salinity_Treatment_Numeric <- as.numeric(as.character(Brach_7day_salinityexposure$Salinity_Treatment)) # Create a dataframe that will store the unique breaks and labels for the x-axis breaks_and_labels <- data.frame( Breaks = sort(unique(Brach_7day_salinityexposure$Salinity_Treatment_Numeric)), Labels = levels(Brach_7day_salinityexposure$Salinity_Treatment_Factor) ) # If there are more breaks than labels, this will truncate the breaks to match the number of labels breaks_and_labels <- breaks_and_labels[1:length(breaks_and_labels$Labels), ] # Plotting Boxplot_Brach_7day_salinityexposure <- ggplot(Brach_7day_salinityexposure, aes(x = Salinity_Treatment_Numeric, y = `Normal (%)`, group = Salinity_Treatment_Factor)) + geom_boxplot(aes(fill = Salinity_Treatment_Factor), show.legend = FALSE, outlier.shape = TRUE, width = 1.5) + scale_x_continuous(name = "Salinity (‰)", breaks = breaks_and_labels$Breaks, labels = breaks_and_labels$Labels) + scale_y_continuous(name = "Normal (%)", limits = c(0, 100)) + scale_fill_manual(values = rep("#bad6f9", length(unique(Brach_4day_salinityexposure$Salinity_Treatment_Factor)))) + ggtitle("7 d") + theme(plot.title = element_text(hjust = 0.5, size = 16, family = "Times New Roman", face = "bold"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), axis.text = element_text(size = 12, family = "Times New Roman"), axis.title = element_text(size = 16, family = "Times New Roman")) # Display the plot print(Boxplot_Brach_7day_salinityexposure) #==================================== # Overlay the model onto the box plot #==================================== Brach_breakpoint_curves_Brach_7day_salinityexposure <- plot_grid(Brach_7day_salinityexposure, breakpoint_curve_Brach_7day_salinityexposure, ncol = 2, nrow = 1) Brach_bps <- grid.arrange(Brach_breakpoint_curves_Brach_7day_salinityexposure) # Ensure that pred_data Salinity_Treatment is numeric and matches the Datas numeric scale pred_data_Brach_7day_salinityexposure$Salinity_Treatment_Numeric <- as.numeric(as.character(pred_data_Brach_7day_salinityexposure$Salinity_Treatment)) Figure_3c <- ggplot(Brach_7day_salinityexposure, aes(x = Salinity_Treatment_Numeric, y = `Normal (%)`)) + geom_boxplot(aes(fill = factor(Salinity_Treatment_Numeric)), show.legend = FALSE, outlier.shape = TRUE, width = 1.5) + scale_fill_manual(values = rep("#f0f0f0", length(unique(Brach_7day_salinityexposure$Salinity_Treatment_Numeric)))) + scale_x_continuous(name = "", breaks = sort(unique(Brach_7day_salinityexposure$Salinity_Treatment_Numeric)), labels = sort(unique(Brach_7day_salinityexposure$Salinity_Treatment_Numeric))) + scale_y_continuous( name = "", limits = c(0, 100), sec.axis = sec_axis(~ . / 100, name = "") # Define secondary axis here ) + geom_line(data = pred_data_Brach_7day_salinityexposure, aes(x = Salinity_Treatment_Numeric, y = Normal_Proportion_Predicted * 100), color = "black", size = 0.5) + geom_vline(xintercept = 22.2, linetype = "dashed", color = "black", size = 0.7) + geom_segment(aes(x = 31.1, y = 5, xend = 31.1, yend = 0), arrow = arrow(type = "closed", length = unit(0.2, "inches")), colour = "#525252", size = 0.5) + # Arrow at 27.9 pointing downwards geom_segment(aes(x = 22.2, y = 5, xend = 22.2, yend = 0), arrow = arrow(type = "closed", length = unit(0.2, "inches")), colour = "#bdbdbd", size = 0.5) + # Arrow at 19 pointing downwards ggtitle("7 days in salinity", "(c)") + theme( plot.title = element_text(hjust = 0.5, size = 30, family = "Times New Roman", face = "bold"), plot.subtitle = element_text(hjust = 0, size = 30, family = "Times New Roman", face = "bold"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), axis.text = element_text(size = 30, family = "Times New Roman"), axis.title = element_text(size = 30, family = "Times New Roman"), axis.title.y.right = element_text(color = "black") ) # Print the plot with arrows pointing downwards print(Figure_3c) greyscale_colors <- c("#f0f0f0", "#d9d9d9", "#bdbdbd", "#969696", "#737373", "#525252", "#252525") ``` #### Survival ```{r} #================================= # Model, pseudo r^2, and LS50/Sopt #================================= #Breakpoint detection Brach_7day_salinityexposure$Surv_prop <- Brach_7day_salinityexposure$Survival / Brach_7day_salinityexposure$`Startin n` model_7day_salinityexposure <- segmented( glm(Surv_prop ~ Salinity_Treatment, data = Brach_7day_salinityexposure, family = binomial(link = "logit")), seg.Z = ~ Salinity_Treatment, psi = 25 ) breakpoint_value <- model_7day_salinityexposure$psi[2] print(breakpoint_value) summary(model_7day_salinityexposure) # pseudo-rsqr logLik_model = logLik(model_7day_salinityexposure) null_model = glm(Surv_prop ~ 1, data = Brach_7day_salinityexposure, family = binomial(link = "logit")) logLik_null = logLik(null_model) R2_McFadden = 1 - (logLik_model / logLik_null) print(R2_McFadden) # LS50 & Sopt90 new_data <- data.frame(Salinity_Treatment = seq(min(Data$Salinity_Treatment), max(Data$Salinity_Treatment), length.out = 1000)) new_data$predicted_prob <- predict(model_7day_salinityexposure, newdata = new_data, type = "response") # Predict probabilities using the segmented model target_probability <- 0.9 closest_index <- which.min(abs(new_data$predicted_prob - target_probability)) # Find the salinity treatment dose where the predicted probability is closest to the target probability dose_at_target_prob <- new_data$Salinity_Treatment[closest_index] print(dose_at_target_prob) target_probability <- 0.5 closest_index <- which.min(abs(new_data$predicted_prob - target_probability)) dose_at_target_prob <- new_data$Salinity_Treatment[closest_index] print(dose_at_target_prob) # =========== # Assumptions # =========== # 1. Linearity of the Logit visreg(model_7day_salinityexposure, "Salinity_Treatment", scale = "response", gg = TRUE) # 2. Goodness-of-Fit - Use Hosmer-Lemeshow test for goodness-of-fit hoslem.test(Brach_7day_salinityexposure$Surv_prop, fitted(model_7day_salinityexposure), g = 10) # 3. Check for Overdispersion overdispersion <- sum(residuals(model_7day_salinityexposure, type = "pearson")^2) / df.residual(model_7day_salinityexposure) print(overdispersion) # 4. Diagnostic plots for residuals par(mfrow = c(2, 2)) plot(model_7day_salinityexposure) ``` ##### Create the prediction curve (breakpoint) and boxplot overlay ```{r} #======================= # Create the model curve #======================= # Generate a sequence of salinity treatments for predictions within the observed range salinity_range_Brach_7day_salinityexposure <- range(Brach_7day_salinityexposure$Salinity_Treatment, na.rm = TRUE) salinity_seq_Brach_7day_salinityexposure <- seq(from = salinity_range_Brach_7day_salinityexposure[1], to = salinity_range_Brach_7day_salinityexposure[2], by = 0.1) # Create a new data frame for predictions pred_data_Brach_7day_salinityexposure <- data.frame(Salinity_Treatment = salinity_seq_Brach_7day_salinityexposure) # Predict survival proportions using the model for 4 day salinity exposure data pred_data_Brach_7day_salinityexposure$Surv_Proportion_Predicted <- predict(model_7day_salinityexposure, newdata = pred_data_Brach_7day_salinityexposure, type = "response") # Create the prediction curve (breakpoint) breakpoint_curve_Brach_7day_salinityexposure <- ggplot(data = pred_data_Brach_7day_salinityexposure) + geom_line(aes(x = Salinity_Treatment, y = Surv_Proportion_Predicted), color = "red", size = 1) + geom_vline(xintercept = 22, linetype = "dashed", color = "grey", size = 1) + scale_x_continuous(breaks = c(17, 19, 22, 25, 30, 34)) + # Custom breaks on numeric x-axis scale_y_continuous(breaks = seq(0, 1, by = 0.1)) + # More tick marks on y-axis xlab("Salinity ‰") + ylab("Prob(Surv)") + ggtitle("7 d") + theme( plot.title = element_text(hjust = 0.5, size = 16, family = "Times New Roman", face = "bold"), panel.grid.major = element_line(color = "grey", size = 0.5), # Add major grid lines panel.grid.minor = element_line(color = "grey", size = 0.25), # Add minor grid lines panel.background = element_rect(fill = "white"), # Set background color to white axis.line = element_line(color = "black"), axis.text = element_text(size = 12, family = "Times New Roman"), axis.title.x = element_text(size = 16, family = "Times New Roman"), axis.title.y = element_text(size = 16, family = "Times New Roman") ) breakpoint_curve_Brach_7day_salinityexposure #========= # Box plot #========= # Assuming 'Salinity_Treatment' is the correct column that needs to be factored and converted to numeric Brach_7day_salinityexposure$Salinity_Treatment_Factor <- factor(Brach_7day_salinityexposure$Salinity_Treatment) Brach_7day_salinityexposure$Salinity_Treatment_Numeric <- as.numeric(as.character(Brach_7day_salinityexposure$Salinity_Treatment)) # Create a dataframe that will store the unique breaks and labels for the x-axis breaks_and_labels <- data.frame( Breaks = sort(unique(Brach_7day_salinityexposure$Salinity_Treatment_Numeric)), Labels = levels(Brach_7day_salinityexposure$Salinity_Treatment_Factor) ) # If there are more breaks than labels, this will truncate the breaks to match the number of labels breaks_and_labels <- breaks_and_labels[1:length(breaks_and_labels$Labels), ] # Plotting Boxplot_Brach_7day_salinityexposure <- ggplot(Brach_7day_salinityexposure, aes(x = Salinity_Treatment_Numeric, y = `Survival (%)`, group = Salinity_Treatment_Factor)) + geom_boxplot(aes(fill = Salinity_Treatment_Factor), show.legend = FALSE, outlier.shape = TRUE, width = 1.5) + scale_x_continuous(name = "Salinity (‰)", breaks = breaks_and_labels$Breaks, labels = breaks_and_labels$Labels) + scale_y_continuous(name = "Survival (%)", limits = c(0, 100)) + scale_fill_manual(values = rep("#bad6f9", length(unique(Brach_4day_salinityexposure$Salinity_Treatment_Factor)))) + ggtitle("7 d") + theme(plot.title = element_text(hjust = 0.5, size = 16, family = "Times New Roman", face = "bold"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), axis.text = element_text(size = 12, family = "Times New Roman"), axis.title = element_text(size = 16, family = "Times New Roman")) # Display the plot print(Boxplot_Brach_7day_salinityexposure) #==================================== # Overlay the model onto the box plot #==================================== Brach_breakpoint_curves_Brach_7day_salinityexposure <- plot_grid(Brach_7day_salinityexposure, breakpoint_curve_Brach_7day_salinityexposure, ncol = 2, nrow = 1) Brach_bps <- grid.arrange(Brach_breakpoint_curves_Brach_7day_salinityexposure) # Ensure that pred_data Salinity_Treatment is numeric and matches the Datas numeric scale pred_data_Brach_7day_salinityexposure$Salinity_Treatment_Numeric <- as.numeric(as.character(pred_data_Brach_7day_salinityexposure$Salinity_Treatment)) Figure_3f <- ggplot(Brach_7day_salinityexposure, aes(x = Salinity_Treatment_Numeric, y = `Survival (%)`)) + geom_boxplot(aes(fill = factor(Salinity_Treatment_Numeric)), show.legend = FALSE, outlier.shape = TRUE, width = 1.5) + scale_fill_manual(values = rep("#f0f0f0", length(unique(Brach_7day_salinityexposure$Salinity_Treatment_Numeric)))) + scale_x_continuous(name = "", breaks = sort(unique(Brach_7day_salinityexposure$Salinity_Treatment_Numeric)), labels = sort(unique(Brach_7day_salinityexposure$Salinity_Treatment_Numeric))) + scale_y_continuous( name = "", limits = c(0, 100), sec.axis = sec_axis(~ . / 100, name = "") ) + geom_line(data = pred_data_Brach_7day_salinityexposure, aes(x = Salinity_Treatment_Numeric, y = Surv_Proportion_Predicted * 100), color = "black", size = 0.5) + geom_vline(xintercept = 22, linetype = "dashed", color = "black", size = 0.7) + geom_segment(aes(x = 28.4, y = 5, xend = 28.4, yend = 0), arrow = arrow(type = "closed", length = unit(0.2, "inches")), colour = "#525252", size = 0.5) + geom_segment(aes(x = 21.9, y = 5, xend = 21.9, yend = 0), arrow = arrow(type = "closed", length = unit(0.2, "inches")), colour = "#bdbdbd", size = 0.5) + ggtitle("", "(f)") + theme( plot.title = element_text(hjust = 0.5, size = 30, family = "Times New Roman", face = "bold"), plot.subtitle = element_text(hjust = 0, size = 30, family = "Times New Roman", face = "bold"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), axis.text = element_text(size = 30, family = "Times New Roman"), axis.title = element_text(size = 30, family = "Times New Roman"), axis.title.y.right = element_text(color = "black") ) # Print the plot with arrows pointing downwards print(Figure_3f) greyscale_colors <- c("#f0f0f0", "#d9d9d9", "#bdbdbd", "#969696", "#737373", "#525252", "#252525") ``` ### Combine the plots ```{r} Figure3 <- plot_grid( Figure_3a, Figure_3b, Figure_3c, Figure_3d, Figure_3e, Figure_3f, ncol = 3, nrow = 2, align = "hv" ) Figure3 ``` # ================================== # Salinity performance curves (SPCs) # ================================== # GLMs ## Data preparation ```{r} Data_stats<-bind_rows( Brach_2day_salinityexposure, Brach_4day_salinityexposure, Brach_7day_salinityexposure ) Data_stats$Sal<-as.factor(Data_stats$Salinity_Treatment) Data_stats$Time<-as.factor(Data_stats$`Pre-settlement assay salinity exposure duration (days)`) Data_stats$cov<-as.factor(Data_stats$`Startin n`) str(Data_stats[, c("Sal", "Time", "cov", "Normal", "Startin n")]) ``` ### Normal larvae model ```{r} # ======== # Table 3a # ======== # response matrix y1 <- cbind(Data_stats$Normal, Data_stats$`Startin n` - Data_stats$Normal) # bias-reduced binomial GLM mod2<-glm(y1 ~ Sal*Time+cov, data = Data_stats, family = binomial(link = "logit"), method = "brglmFit") # Type II test stats::anova(mod2) Sal_Time_emm<-emmeans(mod2, ~Sal*Time, type="response", adjust = "tukey") #USE #Sal_Time_emm<-emmeans(mod2, ~Sal, type="response", adjust = "tukey") #USE pairs <- pairs(Sal_Time_emm) pairs # ===================== # Supplementary Table 3 # ===================== # Convert the pairs output to a data frame if necessary pairs_df <- as.data.frame(pairs) normtable <- kable(pairs_df, format = "html", table.attr = 'border="0"') writeLines(normtable, '[insert filepath]/Supp table 3.xls') # ===================== # Figure 4a # ===================== emmip(mod2, ~Sal|Time, type="response", CI=TRUE, adjust = "tukey") #USE emmip(mod2, ~Sal, type="response", CI=TRUE, adjust = "tukey") #test posthocdf2<-as.data.frame(Sal_Time_emm[1:18]) posthocdf2[1,7]<-0 posthocdf2[7,7]<-0 posthocdf2[13,7]<-0 posthocdf2[14,7]<-0 Figure_4a <- ggplot(posthocdf2, aes(x = Sal, y = prob, group = Time, color = Time)) + geom_line(position = position_dodge(width = 0.25)) + geom_point(position = position_dodge(width = 0.25)) + geom_ribbon(aes(ymin = asymp.LCL, ymax = asymp.UCL, fill = Time), alpha = 0.2, position = position_dodge(width = 0.25)) + xlab("Salinity (‰)") + ylab("Prob(Normal)") + ggtitle("", "(a)") + theme( plot.title = element_text(hjust = 0.5, size = 16, family = "Times New Roman", face = "bold"), plot.subtitle = element_text(size = 16, hjust = -0.05, vjust = 5, family = "Times New Roman", face = "bold"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), axis.text = element_text(size = 16, family = "Times New Roman"), strip.background = element_blank(), strip.text = element_text(size = 16, family = "Times New Roman"), axis.title.y = element_text(size = 16, family = "Times New Roman"), axis.title.x = element_text(size = 16, family = "Times New Roman"), legend.position = "", legend.text = element_text(size = 16, family = "Times New Roman"), legend.title = element_text(size = 16, family = "Times New Roman") ) + labs(color = "Larval salinity treatment \nduration (days)", fill = "Larval salinity treatment \nduration (days)") + scale_color_grey(start = 0.85, end = 0) + scale_fill_grey(start = 0.85, end = 0) print(Figure_4a) summary(mod2) exp(coef(mod2)) mod2$fitted.values tidy(mod2) conf_intervals <- confint(mod2, level = 0.95) print(conf_intervals) ``` ### Survival model ```{r} # ======== # Table 3b # ======== # response matrix y2 <- cbind(Data_stats$Survival, Data_stats$`Startin n` - Data_stats$Survival) # bias-reduced binomial GLM mod3<-glm(y2 ~ Sal*Time+cov, data = Data_stats, family = binomial(link = "logit"), method = "brglmFit") # Type II test stats::anova(mod3) Sal_Time_emm<-emmeans(mod3, ~Sal*Time, type="response", adjust = "tukey") #USE #Sal_Time_emm<-emmeans(mod3, ~Sal, type="response", adjust = "tukey") #USE pairs <- pairs(Sal_Time_emm) pairs # ===================== # Supplementary Table 4 # ===================== # Convert the pairs output to a data frame if necessary pairs_df <- as.data.frame(pairs) normtable <- kable(pairs_df, format = "html", table.attr = 'border="0"') writeLines(normtable, '[insert filepath]/Supp table 4.xls') # ===================== # Figure 4b # ===================== emmip(mod3, ~Sal|Time, type="response", CI=TRUE, adjust = "tukey") #USE emmip(mod3, ~Sal, type="response", CI=TRUE, adjust = "tukey") #test posthocdf3<-as.data.frame(Sal_Time_emm[1:18]) posthocdf3[1,7]<-0 posthocdf3[7,7]<-0 posthocdf3[13,7]<-0 posthocdf3[14,7]<-0 Figure_4b <- ggplot(posthocdf3, aes(x = Sal, y = prob, group = Time, color = Time)) + geom_line(position = position_dodge(width = 0.25)) + geom_point(position = position_dodge(width = 0.25)) + geom_ribbon(aes(ymin = asymp.LCL, ymax = asymp.UCL, fill = Time), alpha = 0.2, position = position_dodge(width = 0.25)) + xlab("Salinity (‰)") + ylab("Prob(Survival)") + ggtitle("", "(b)") + theme( plot.title = element_text(hjust = 0.5, size = 16, family = "Times New Roman", face = "bold"), plot.subtitle = element_text(size = 16, hjust = -0.05, vjust = 5, family = "Times New Roman", face = "bold"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), axis.text = element_text(size = 16, family = "Times New Roman"), strip.background = element_blank(), strip.text = element_text(size = 16, family = "Times New Roman"), axis.title.y = element_text(size = 16, family = "Times New Roman"), axis.title.x = element_text(size = 16, family = "Times New Roman"), legend.position = "", legend.text = element_text(size = 16, family = "Times New Roman"), legend.title = element_text(size = 16, family = "Times New Roman") ) + labs(color = "Larval salinity treatment \nduration (days)", fill = "Larval salinity treatment \nduration (days)") + scale_color_grey(start = 0.85, end = 0) + scale_fill_grey(start = 0.85, end = 0) print(Figure_4b) summary(mod3) exp(coef(mod3)) mod2$fitted.values tidy(mod3) conf_intervals <- confint(mod3, level = 0.95) print(conf_intervals) ``` ### Combine the plots ```{r} Figure4 <- plot_grid( Figure_4a, Figure_4b, ncol = 2, align = "hv" ) Figure4 ``` # ======================== # Larval settlement assays # ======================== ### Metamorphosis and juvenile success mean summary #### Load data ```{r} Set_Data_long<-read_excel('[insert filepath]/Clements & Byrne 2026.xlsx', sheet = "Larval settlement assays (long)", range = "A2:E259", col_names = c("Low salinity exposure duration (d)", "Duration in settlement assay", "Salinity", "Metamorphic response", "Value (%)")) Set_Data$Salinity<-as.factor(Set_Data$`Larval Salinity Group`) Set_Data$Time_exposed<-as.factor(Set_Data$`Pre-settlement assay salinity exposure duration (days)`) Set_Data$Time_in_assay<-as.factor(Set_Data$`Duration in settlement assay (days)`) ``` #### 2 day larval salinity exposure ```{r} # ===================== # Supplementary table 5 # ===================== Set_Data_2d_2d <- Set_Data %>% filter(Time_exposed == "2", Time_in_assay == "2") Set_2d_2d <- Set_Data_2d_2d %>% filter(`Larval Salinity Group` %in% c(25, 30, 34)) %>% group_by(`Larval Salinity Group`) %>% summarise( Mean = mean(`Settlement (%)`, na.rm = TRUE), Std_Error = sd(`Settlement (%)`, na.rm = TRUE) / sqrt(n()), count = n()) Set_2d_2d Set_2d_2d_mm <- Set_Data_2d_2d %>% filter(`Larval Salinity Group` %in% c(25, 30, 34)) %>% group_by(`Larval Salinity Group`) %>% summarise( Mean = mean(`Undergoing metamorphosis (%)`, na.rm = TRUE), Std_Error = sd(`Undergoing metamorphosis (%)`, na.rm = TRUE) / sqrt(n()), count = n()) Set_2d_2d_mm Set_2d_2d_J <- Set_Data_2d_2d %>% filter(`Larval Salinity Group` %in% c(25, 30, 34)) %>% group_by(`Larval Salinity Group`) %>% summarise( Mean = mean(`Juveniles (%)`, na.rm = TRUE), Std_Error = sd(`Juveniles (%)`, na.rm = TRUE) / sqrt(n()), count = n()) Set_2d_2d_J Set_Data_2d_5d <- Set_Data %>% filter(Time_exposed == "2", Time_in_assay == "5") Set_2d_5d <- Set_Data_2d_5d %>% filter(`Larval Salinity Group` %in% c(25, 30, 34)) %>% group_by(`Larval Salinity Group`) %>% summarise( Mean = mean(`Settlement (%)`, na.rm = TRUE), Std_Error = sd(`Settlement (%)`, na.rm = TRUE) / sqrt(n()), count = n()) Set_2d_5d Set_2d_5d_mm <- Set_Data_2d_5d %>% filter(`Larval Salinity Group` %in% c(25, 30, 34)) %>% group_by(`Larval Salinity Group`) %>% summarise( Mean = mean(`Undergoing metamorphosis (%)`, na.rm = TRUE), Std_Error = sd(`Undergoing metamorphosis (%)`, na.rm = TRUE) / sqrt(n()), count = n()) Set_2d_5d_mm Set_2d_5d_J <- Set_Data_2d_5d %>% filter(`Larval Salinity Group` %in% c(25, 30, 34)) %>% group_by(`Larval Salinity Group`) %>% summarise( Mean = mean(`Juveniles (%)`, na.rm = TRUE), Std_Error = sd(`Juveniles (%)`, na.rm = TRUE) / sqrt(n()), count = n()) Set_2d_5d_J Set_Data_2d_7d <- Set_Data %>% filter(Time_exposed == "2", Time_in_assay == "7") Set_2d_7d <- Set_Data_2d_7d %>% filter(`Larval Salinity Group` %in% c(25, 30, 34)) %>% group_by(`Larval Salinity Group`) %>% summarise( Mean = mean(`Settlement (%)`, na.rm = TRUE), Std_Error = sd(`Settlement (%)`, na.rm = TRUE) / sqrt(n()), count = n()) Set_2d_7d Set_2d_7d_mm <- Set_Data_2d_7d %>% filter(`Larval Salinity Group` %in% c(25, 30, 34)) %>% group_by(`Larval Salinity Group`) %>% summarise( Mean = mean(`Undergoing metamorphosis (%)`, na.rm = TRUE), Std_Error = sd(`Undergoing metamorphosis (%)`, na.rm = TRUE) / sqrt(n()), count = n()) Set_2d_7d_mm Set_2d_7d_J <- Set_Data_2d_7d %>% filter(`Larval Salinity Group` %in% c(25, 30, 34)) %>% group_by(`Larval Salinity Group`) %>% summarise( Mean = mean(`Juveniles (%)`, na.rm = TRUE), Std_Error = sd(`Juveniles (%)`, na.rm = TRUE) / sqrt(n()), count = n()) Set_2d_7d_J ``` #### 4 day larval salinity exposure ```{r} # ===================== # Supplementary table 5 # ===================== #split & reshape the data Set_Data_4d_2d <- Set_Data %>% filter(Time_exposed == "4", Time_in_assay == "2") Set_4d_2d <- Set_Data_4d_2d %>% filter(`Larval Salinity Group` %in% c(25, 30, 34)) %>% group_by(`Larval Salinity Group`) %>% summarise( Mean = mean(`Settlement (%)`, na.rm = TRUE), Std_Error = sd(`Settlement (%)`, na.rm = TRUE) / sqrt(n()), count = n()) Set_4d_2d Set_4d_2d_mm <- Set_Data_4d_2d %>% filter(`Larval Salinity Group` %in% c(25, 30, 34)) %>% group_by(`Larval Salinity Group`) %>% summarise( Mean = mean(`Undergoing metamorphosis (%)`, na.rm = TRUE), Std_Error = sd(`Undergoing metamorphosis (%)`, na.rm = TRUE) / sqrt(n()), count = n()) Set_4d_2d_mm Set_4d_2d_J <- Set_Data_4d_2d %>% filter(`Larval Salinity Group` %in% c(25, 30, 34)) %>% group_by(`Larval Salinity Group`) %>% summarise( Mean = mean(`Juveniles (%)`, na.rm = TRUE), Std_Error = sd(`Juveniles (%)`, na.rm = TRUE) / sqrt(n()), count = n()) Set_4d_2d_J Set_Data_4d_5d <- Set_Data %>% filter(Time_exposed == "4", Time_in_assay == "5") Set_4d_5d <- Set_Data_4d_5d %>% filter(`Larval Salinity Group` %in% c(25, 30, 34)) %>% group_by(`Larval Salinity Group`) %>% summarise( Mean = mean(`Settlement (%)`, na.rm = TRUE), Std_Error = sd(`Settlement (%)`, na.rm = TRUE) / sqrt(n()), count = n()) Set_4d_5d Set_4d_5d_mm <- Set_Data_4d_5d %>% filter(`Larval Salinity Group` %in% c(25, 30, 34)) %>% group_by(`Larval Salinity Group`) %>% summarise( Mean = mean(`Undergoing metamorphosis (%)`, na.rm = TRUE), Std_Error = sd(`Undergoing metamorphosis (%)`, na.rm = TRUE) / sqrt(n()), count = n()) Set_4d_5d_mm Set_4d_5d_J <- Set_Data_4d_5d %>% filter(`Larval Salinity Group` %in% c(25, 30, 34)) %>% group_by(`Larval Salinity Group`) %>% summarise( Mean = mean(`Juveniles (%)`, na.rm = TRUE), Std_Error = sd(`Juveniles (%)`, na.rm = TRUE) / sqrt(n()), count = n()) Set_4d_5d_J Set_Data_4d_7d <- Set_Data %>% filter(Time_exposed == "4", Time_in_assay == "7") Set_4d_7d <- Set_Data_4d_7d %>% filter(`Larval Salinity Group` %in% c(25, 30, 34)) %>% group_by(`Larval Salinity Group`) %>% summarise( Mean = mean(`Settlement (%)`, na.rm = TRUE), Std_Error = sd(`Settlement (%)`, na.rm = TRUE) / sqrt(n()), count = n()) Set_4d_7d Set_4d_7d_mm <- Set_Data_4d_7d %>% filter(`Larval Salinity Group` %in% c(25, 30, 34)) %>% group_by(`Larval Salinity Group`) %>% summarise( Mean = mean(`Undergoing metamorphosis (%)`, na.rm = TRUE), Std_Error = sd(`Undergoing metamorphosis (%)`, na.rm = TRUE) / sqrt(n()), count = n()) Set_4d_7d_mm Set_4d_7d_J <- Set_Data_4d_7d %>% filter(`Larval Salinity Group` %in% c(25, 30, 34)) %>% group_by(`Larval Salinity Group`) %>% summarise( Mean = mean(`Juveniles (%)`, na.rm = TRUE), Std_Error = sd(`Juveniles (%)`, na.rm = TRUE) / sqrt(n()), count = n()) Set_4d_7d_J ``` #### 7 day larval salinity exposure ```{r} # ===================== # Supplementary table 5 # ===================== #split & reshape the data Set_Data_7d_2d <- Set_Data %>% filter(Time_exposed == "7", Time_in_assay == "2") Set_7d_2d <- Set_Data_7d_2d %>% filter(`Larval Salinity Group` %in% c(25, 30, 34)) %>% group_by(`Larval Salinity Group`) %>% summarise( Mean = mean(`Settlement (%)`, na.rm = TRUE), Std_Error = sd(`Settlement (%)`, na.rm = TRUE) / sqrt(n()), count = n()) Set_7d_2d Set_Data_7d_2d_mm <- Set_Data_7d_2d %>% filter(`Larval Salinity Group` %in% c(25, 30, 34)) %>% group_by(`Larval Salinity Group`) %>% summarise( Mean = mean(`Undergoing metamorphosis (%)`, na.rm = TRUE), Std_Error = sd(`Undergoing metamorphosis (%)`, na.rm = TRUE) / sqrt(n()), count = n()) Set_Data_7d_2d_mm Set_7d_2d_J <- Set_Data_7d_2d %>% filter(`Larval Salinity Group` %in% c(25, 30, 34)) %>% group_by(`Larval Salinity Group`) %>% summarise( Mean = mean(`Juveniles (%)`, na.rm = TRUE), Std_Error = sd(`Juveniles (%)`, na.rm = TRUE) / sqrt(n()), count = n()) Set_7d_2d_J Set_Data_7d_5d <- Set_Data %>% filter(Time_exposed == "7", Time_in_assay == "5") Set_7d_5d <- Set_Data_7d_5d %>% filter(`Larval Salinity Group` %in% c(25, 30, 34)) %>% group_by(`Larval Salinity Group`) %>% summarise( Mean = mean(`Settlement (%)`, na.rm = TRUE), Std_Error = sd(`Settlement (%)`, na.rm = TRUE) / sqrt(n()), count = n()) Set_7d_5d Set_Data_7d_5d_mm <- Set_Data_7d_5d %>% filter(`Larval Salinity Group` %in% c(25, 30, 34)) %>% group_by(`Larval Salinity Group`) %>% summarise( Mean = mean(`Undergoing metamorphosis (%)`, na.rm = TRUE), Std_Error = sd(`Undergoing metamorphosis (%)`, na.rm = TRUE) / sqrt(n()), count = n()) Set_Data_7d_5d_mm Set_7d_5d_J <- Set_Data_7d_5d %>% filter(`Larval Salinity Group` %in% c(25, 30, 34)) %>% group_by(`Larval Salinity Group`) %>% summarise( Mean = mean(`Juveniles (%)`, na.rm = TRUE), Std_Error = sd(`Juveniles (%)`, na.rm = TRUE) / sqrt(n()), count = n()) Set_7d_5d_J Set_Data_7d_7d <- Set_Data %>% filter(Time_exposed == "7", Time_in_assay == "7") Set_7d_7d <- Set_Data_7d_7d %>% filter(`Larval Salinity Group` %in% c(25, 30, 34)) %>% group_by(`Larval Salinity Group`) %>% summarise( Mean = mean(`Settlement (%)`, na.rm = TRUE), Std_Error = sd(`Settlement (%)`, na.rm = TRUE) / sqrt(n()), count = n()) Set_7d_7d Set_Data_7d_7d_mm <- Set_Data_7d_7d %>% filter(`Larval Salinity Group` %in% c(25, 30, 34)) %>% group_by(`Larval Salinity Group`) %>% summarise( Mean = mean(`Undergoing metamorphosis (%)`, na.rm = TRUE), Std_Error = sd(`Undergoing metamorphosis (%)`, na.rm = TRUE) / sqrt(n()), count = n()) Set_Data_7d_7d_mm Set_7d_7d_J <- Set_Data_7d_7d %>% filter(`Larval Salinity Group` %in% c(25, 30, 34)) %>% group_by(`Larval Salinity Group`) %>% summarise( Mean = mean(`Juveniles (%)`, na.rm = TRUE), Std_Error = sd(`Juveniles (%)`, na.rm = TRUE) / sqrt(n()), count = n()) Set_7d_7d_J ``` ##### Stacked Column graphs ###### Format data for analyses ```{r} Set_Data_long$Salinity<-as.factor(Set_Data_long$Salinity) Set_Data_long$Time_exposed<-as.factor(Set_Data_long$`Low salinity exposure duration (d)`) Set_Data_long$Time_in_assay<-as.factor(Set_Data_long$`Duration in settlement assay`) y<-Set_Data_long$`Value (%)` ``` ###### 2d sal tol metamorphic response ```{r} # ======== # Figure 5 # ======== # Create the ggplot with cumulative SE greyscale_colors <- c("#f0f0f0", "#d9d9d9", "#bdbdbd", "#969696", "#737373", "#525252", "#252525") #2 d Set_Data_long_2d_2d<- Set_Data_long %>% filter(`Low salinity exposure duration (d)`=="2",`Duration in settlement assay` == "2") Set_Data_long_2d_2d$`Metamorphic response`<-as.factor(Set_Data_long_2d_2d$`Metamorphic response`) df3 <- Set_Data_long_2d_2d %>% group_by(Salinity, `Metamorphic response`) %>% summarise(mean3 = mean(`Value (%)`), se3 = sd(`Value (%)`) / sqrt(n()), # Calculate standard error .groups = 'drop') %>% mutate(Salinity = factor(Salinity, levels = c(25, 30, 34)), `Metamorphic response` = factor(`Metamorphic response`)) df3$`Metamorphic response` <- factor(df3$`Metamorphic response`, levels = c("Metamorphosis", "Juveniles")) # Calculate the position for error bars error_bars3 = df3 %>% arrange(Salinity, desc(`Metamorphic response`)) %>% group_by(Salinity) %>% mutate(mean3 = lag(cumsum(mean3), default = 0) + mean3) %>% ungroup() p2 <- ggplot(df3, aes(x = Salinity, y = mean3, fill = `Metamorphic response`)) + geom_bar(stat = 'identity') + geom_errorbar(data = error_bars3, aes(ymax = mean3 + se3, ymin = mean3 - se3), position = position_dodge(width = 0.3), width = 0.2) + ggtitle("2d low salinity exposure","2 days in settlement assay") + labs(x = "", y = "", fill = NULL) + theme(plot.title = element_text(hjust = 0.0, vjust = 1, size = 24, family = "Times New Roman", face = "bold"), plot.subtitle = element_text(hjust = 0.5, size = 24, family = "Times New Roman"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), axis.text = element_text(size = 30, family = "Times New Roman"), axis.title.y = element_text(size = 30, family = "Times New Roman"), axis.title.x = element_text(size = 30, family = "Times New Roman"), legend.position = "") + # Adjust legend position to the right scale_fill_manual(values = c("#d9d9d9","#737373")) + ylim(0, 120) + scale_y_continuous(limits = c(0, 120),breaks = c(0, 20, 40, 60, 80, 100)) p2 Set_Data_long_2d_5d<- Set_Data_long %>% filter(`Low salinity exposure duration (d)`=="2",`Duration in settlement assay` == "5") Set_Data_long_2d_5d$`Metamorphic response`<-as.factor(Set_Data_long_2d_5d$`Metamorphic response`) df6 <- Set_Data_long_2d_5d %>% group_by(Salinity, `Metamorphic response`) %>% summarise(mean6 = mean(`Value (%)`), se6 = sd(`Value (%)`) / sqrt(n()), # Calculate standard error .groups = 'drop') %>% mutate(Salinity = factor(Salinity, levels = c(25, 30, 34)), `Metamorphic response` = factor(`Metamorphic response`)) df6$`Metamorphic response` <- factor(df6$`Metamorphic response`, levels = c("Metamorphosis", "Juveniles")) # Calculate the position for error bars error_bars6 = df6 %>% arrange(Salinity, desc(`Metamorphic response`)) %>% group_by(Salinity) %>% mutate(mean6 = lag(cumsum(mean6), default = 0) + mean6) %>% ungroup() p3 <- ggplot(df6, aes(x = Salinity, y = mean6, fill = `Metamorphic response`)) + geom_bar(stat = 'identity') + geom_errorbar(data = error_bars6, aes(ymax = mean6 + se6, ymin = mean6 - se6), position = position_dodge(width = 0.3), width = 0.2) + ggtitle("","5 days in settlement assay") + labs(x = "", y = "", fill = NULL) + theme(plot.title = element_text(hjust = 0, size = 24, family = "Times New Roman", face = "bold"), plot.subtitle = element_text(hjust = 0.5, size = 24, family = "Times New Roman"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), axis.text = element_text(size = 30, family = "Times New Roman"), axis.title.y = element_text(size = 30, family = "Times New Roman"), axis.title.x = element_text(size = 30, family = "Times New Roman"), legend.position = "") + # Adjust legend position to the right scale_fill_manual(values = c("#d9d9d9","#737373")) + ylim(0, 120) + scale_y_continuous(limits = c(0, 120),breaks = c(0, 20, 40, 60, 80, 100)) # Set y-axis limits p3 Set_Data_long_2d_7d<- Set_Data_long %>% filter(`Low salinity exposure duration (d)`=="2",`Duration in settlement assay` == "7") Set_Data_long_2d_7d$`Metamorphic response`<-as.factor(Set_Data_long_2d_7d$`Metamorphic response`) df8 <- Set_Data_long_2d_7d %>% group_by(Salinity, `Metamorphic response`) %>% summarise(mean8 = mean(`Value (%)`), se8 = sd(`Value (%)`) / sqrt(n()), # Calculate standard error .groups = 'drop') %>% mutate(Salinity = factor(Salinity, levels = c(25, 30, 34)), `Metamorphic response` = factor(`Metamorphic response`)) df8$`Metamorphic response` <- factor(df8$`Metamorphic response`, levels = c("Metamorphosis", "Juveniles")) # Calculate the position for error bars error_bars8 = df8 %>% arrange(Salinity, desc(`Metamorphic response`)) %>% group_by(Salinity) %>% mutate(mean8 = lag(cumsum(mean8), default = 0) + mean8) %>% ungroup() p4 <- ggplot(df8, aes(x = Salinity, y = mean8, fill = `Metamorphic response`)) + geom_bar(stat = 'identity') + geom_errorbar(data = error_bars8, aes(ymax = mean8 + se8, ymin = mean8 - se8), position = position_dodge(width = 0.3), width = 0.2) + ggtitle("","7 days in settlement assay") + labs(x = "", y = "", fill = NULL) + theme(plot.title = element_text(hjust = 0, size = 24, family = "Times New Roman", face = "bold"), plot.subtitle = element_text(hjust = 0.5, size = 24, family = "Times New Roman"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), axis.text = element_text(size = 30, family = "Times New Roman"), axis.title.y = element_text(size = 30, family = "Times New Roman"), axis.title.x = element_text(size = 30, family = "Times New Roman"), legend.position = "") + # Adjust legend position to the right scale_fill_manual(values = c("#d9d9d9","#737373")) + ylim(0, 120) + scale_y_continuous(limits = c(0, 120),breaks = c(0, 20, 40, 60, 80, 100)) # Set y-axis limits p4 stackedbps_2d<-plot_grid(p2, p3, p4, ncol = 3,nrow = 1) stackedbps_2d<-grid.arrange(arrangeGrob(stackedbps_2d)) ``` ###### 4d sal tol metamorphic response ```{r} # ======== # Figure 5 # ======== #4 d Set_Data_long_4d_2d<- Set_Data_long %>% filter(`Low salinity exposure duration (d)`=="4",`Duration in settlement assay` == "2") Set_Data_long_4d_2d$`Metamorphic response`<-as.factor(Set_Data_long_4d_2d$`Metamorphic response`) df11 <- Set_Data_long_4d_2d %>% group_by(Salinity, `Metamorphic response`) %>% summarise(mean11 = mean(`Value (%)`), se11 = sd(`Value (%)`) / sqrt(n()), # Calculate standard error .groups = 'drop') %>% mutate(Salinity = factor(Salinity, levels = c(25, 30, 34)), `Metamorphic response` = factor(`Metamorphic response`)) df11$`Metamorphic response` <- factor(df11$`Metamorphic response`, levels = c("Metamorphosis", "Juveniles")) # Calculate the position for error bars error_bars11 = df11 %>% arrange(Salinity, desc(`Metamorphic response`)) %>% group_by(Salinity) %>% mutate(mean11 = lag(cumsum(mean11), default = 0) + mean11) %>% ungroup() p10 <- ggplot(df11, aes(x = Salinity, y = mean11, fill = `Metamorphic response`)) + geom_bar(stat = 'identity') + geom_errorbar(data = error_bars11, aes(ymax = mean11 + se11, ymin = mean11 - se11), position = position_dodge(width = 0.3), width = 0.2) + ggtitle("4d low salinity exposure","") + labs(x = "", y = "Metamorphic response (%)", fill = NULL) + theme(plot.title = element_text(hjust = 0, size = 24, family = "Times New Roman", face = "bold"), plot.subtitle = element_text(hjust = 0.5, size = 24, family = "Times New Roman"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), axis.text = element_text(size = 30, family = "Times New Roman"), axis.title.y = element_text(size = 30, family = "Times New Roman"), axis.title.x = element_text(size = 30, family = "Times New Roman"), legend.position = "") + # Adjust legend position to the right scale_fill_manual(values = c("#d9d9d9","#737373")) + ylim(0, 120) + scale_y_continuous(limits = c(0, 120),breaks = c(0, 20, 40, 60, 80, 100)) # Set y-axis limits p10 Set_Data_long_4d_5d<- Set_Data_long %>% filter(`Low salinity exposure duration (d)`=="4",`Duration in settlement assay` == "5") Set_Data_long_4d_5d$`Metamorphic response`<-as.factor(Set_Data_long_4d_5d$`Metamorphic response`) df14 <- Set_Data_long_4d_5d %>% group_by(Salinity, `Metamorphic response`) %>% summarise(mean14 = mean(`Value (%)`), se14 = sd(`Value (%)`) / sqrt(n()), # Calculate standard error .groups = 'drop') %>% mutate(Salinity = factor(Salinity, levels = c(25, 30, 34)), `Metamorphic response` = factor(`Metamorphic response`)) df14$`Metamorphic response` <- factor(df14$`Metamorphic response`, levels = c("Metamorphosis", "Juveniles")) # Calculate the position for error bars error_bars14 = df14 %>% arrange(Salinity, desc(`Metamorphic response`)) %>% group_by(Salinity) %>% mutate(mean14 = lag(cumsum(mean14), default = 0) + mean14) %>% ungroup() p11 <- ggplot(df14, aes(x = Salinity, y = mean14, fill = `Metamorphic response`)) + geom_bar(stat = 'identity') + geom_errorbar(data = error_bars14, aes(ymax = mean14 + se14, ymin = mean14 - se14), position = position_dodge(width = 0.3), width = 0.2) + ggtitle("","") + labs(x = "", y = "", fill = NULL) + theme(plot.title = element_text(hjust = 0, size = 24, family = "Times New Roman", face = "bold"), plot.subtitle = element_text(hjust = 0.5, size = 24, family = "Times New Roman"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), axis.text = element_text(size = 30, family = "Times New Roman"), axis.title.y = element_text(size = 30, family = "Times New Roman"), axis.title.x = element_text(size = 30, family = "Times New Roman"), legend.position = "") + # Adjust legend position to the right scale_fill_manual(values = c("#d9d9d9","#737373")) + ylim(0, 120) + scale_y_continuous(limits = c(0, 120),breaks = c(0, 20, 40, 60, 80, 100)) # Set y-axis limits p11 Set_Data_long_4d_7d<- Set_Data_long %>% filter(`Low salinity exposure duration (d)`=="4",`Duration in settlement assay` == "7") Set_Data_long_4d_7d$`Metamorphic response`<-as.factor(Set_Data_long_4d_7d$`Metamorphic response`) df16 <- Set_Data_long_4d_7d %>% group_by(Salinity, `Metamorphic response`) %>% summarise(mean16 = mean(`Value (%)`), se16 = sd(`Value (%)`) / sqrt(n()), # Calculate standard error .groups = 'drop') %>% mutate(Salinity = factor(Salinity, levels = c(25, 30, 34)), `Metamorphic response` = factor(`Metamorphic response`)) df16$`Metamorphic response` <- factor(df16$`Metamorphic response`, levels = c("Metamorphosis", "Juveniles")) # Calculate the position for error bars error_bars16 = df16 %>% arrange(Salinity, desc(`Metamorphic response`)) %>% group_by(Salinity) %>% mutate(mean16 = lag(cumsum(mean16), default = 0) + mean16) %>% ungroup() p12 <- ggplot(df16, aes(x = Salinity, y = mean16, fill = `Metamorphic response`)) + geom_bar(stat = 'identity') + geom_errorbar(data = error_bars16, aes(ymax = mean16 + se16, ymin = mean16 - se16), position = position_dodge(width = 0.3), width = 0.2) + ggtitle("","") + labs(x = "", y = "", fill = NULL) + theme(plot.title = element_text(hjust = 0, size = 24, family = "Times New Roman", face = "bold"), plot.subtitle = element_text(hjust = 0.5, size = 24, family = "Times New Roman"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), axis.text = element_text(size = 30, family = "Times New Roman"), axis.title.y = element_text(size = 30, family = "Times New Roman"), axis.title.x = element_text(size = 30, family = "Times New Roman"), legend.position = "") + # Adjust legend position to the right scale_fill_manual(values = c("#d9d9d9","#737373")) + ylim(0, 120) + scale_y_continuous(limits = c(0, 120),breaks = c(0, 20, 40, 60, 80, 100)) # Set y-axis limits p12 stackedbps_4d<-plot_grid(p10, p11, p12, ncol = 3,nrow = 1) stackedbps_4d<-grid.arrange(arrangeGrob(stackedbps_4d)) ``` ###### 7d sal tol metamorphic response ```{r} # ======== # Figure 5 # ======== #7 d Set_Data_long_7d_0d<- Set_Data_long %>% filter(`Low salinity exposure duration (d)`=="7",`Duration in settlement assay` == "0") Set_Data_long_7d_0d$`Metamorphic response`<-as.factor(Set_Data_long_7d_0d$`Metamorphic response`) df17 <- Set_Data_long_7d_0d %>% group_by(Salinity, `Metamorphic response`) %>% summarise(mean17 = mean(`Value (%)`), se17 = sd(`Value (%)`) / sqrt(n()), # Calculate standard error .groups = 'drop') %>% mutate(Salinity = factor(Salinity, levels = c(34, 30, 25)), `Metamorphic response` = factor(`Metamorphic response`)) df17$`Metamorphic response` <- factor(df17$`Metamorphic response`, levels = c("Metamorphosis", "Juveniles")) # Calculate the position for error bars error_bars17 = df17 %>% arrange(Salinity, desc(`Metamorphic response`)) %>% group_by(Salinity) %>% mutate(mean17 = lag(cumsum(mean17), default = 0) + mean17) %>% ungroup() p17 <- ggplot(df17, aes(x = Salinity, y = mean17, fill = `Metamorphic response`)) + geom_bar(stat = 'identity') + geom_errorbar(data = error_bars17, aes(ymax = mean17 + se17, ymin = mean17 - se17), position = position_dodge(width = 0.3), width = 0.2) + ggtitle("7d low salinity exposure","") + labs(x = "", y = "", fill = NULL) + theme(plot.title = element_text(hjust = 0, size = 24, family = "Times New Roman", face = "bold"), plot.subtitle = element_text(hjust = 0.5, size = 24, family = "Times New Roman"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), axis.text = element_text(size = 30, family = "Times New Roman"), axis.title.y = element_text(size = 30, family = "Times New Roman"), axis.title.x = element_text(size = 30, family = "Times New Roman"), legend.position = "") + # Adjust legend position to the right scale_fill_manual(values = c("#d9d9d9","#737373")) + ylim(0, 120) + scale_y_continuous(limits = c(0, 120),breaks = c(0, 20, 40, 60, 80, 100)) # Set y-axis limits p17 Set_Data_long_7d_2d<- Set_Data_long %>% filter(`Low salinity exposure duration (d)`=="7",`Duration in settlement assay` == "2") Set_Data_long_7d_2d$`Metamorphic response`<-as.factor(Set_Data_long_7d_2d$`Metamorphic response`) df18 <- Set_Data_long_7d_2d %>% group_by(Salinity, `Metamorphic response`) %>% summarise(mean18 = mean(`Value (%)`), se18 = sd(`Value (%)`) / sqrt(n()), # Calculate standard error .groups = 'drop') %>% mutate(Salinity = factor(Salinity, levels = c(25, 30, 34)), `Metamorphic response` = factor(`Metamorphic response`)) df18$`Metamorphic response` <- factor(df18$`Metamorphic response`, levels = c("Metamorphosis", "Juveniles")) # Calculate the position for error bars error_bars18 = df18 %>% arrange(Salinity, desc(`Metamorphic response`)) %>% group_by(Salinity) %>% mutate(mean18 = lag(cumsum(mean18), default = 0) + mean18) %>% ungroup() p18 <- ggplot(df18, aes(x = Salinity, y = mean18, fill = `Metamorphic response`)) + geom_bar(stat = 'identity') + geom_errorbar(data = error_bars18, aes(ymax = mean18 + se18, ymin = mean18 - se18), position = position_dodge(width = 0.3), width = 0.2) + ggtitle("7d low salinity exposure","") + labs(x = "", y = "", fill = NULL) + theme(plot.title = element_text(hjust = 0, size = 24, family = "Times New Roman", face = "bold"), plot.subtitle = element_text(hjust = 0.5, size = 24, family = "Times New Roman"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), axis.text = element_text(size = 30, family = "Times New Roman"), axis.title.y = element_text(size = 30, family = "Times New Roman"), axis.title.x = element_text(size = 30, family = "Times New Roman"), legend.position = "") + # Adjust legend position to the right scale_fill_manual(values = c("#d9d9d9","#737373")) + ylim(0, 120) + scale_y_continuous(limits = c(0, 120),breaks = c(0, 20, 40, 60, 80, 100)) # Set y-axis limits p18 Set_Data_long_7d_5d<- Set_Data_long %>% filter(`Low salinity exposure duration (d)`=="7",`Duration in settlement assay` == "5") Set_Data_long_7d_5d$`Metamorphic response`<-as.factor(Set_Data_long_7d_5d$`Metamorphic response`) df19 <- Set_Data_long_7d_5d %>% group_by(Salinity, `Metamorphic response`) %>% summarise(mean19 = mean(`Value (%)`), se19 = sd(`Value (%)`) / sqrt(n()), # Calculate standard error .groups = 'drop') %>% mutate(Salinity = factor(Salinity, levels = c(25, 30, 34)), `Metamorphic response` = factor(`Metamorphic response`)) df19$`Metamorphic response` <- factor(df19$`Metamorphic response`, levels = c("Metamorphosis", "Juveniles")) # Calculate the position for error bars error_bars19 = df19 %>% arrange(Salinity, desc(`Metamorphic response`)) %>% group_by(Salinity) %>% mutate(mean19 = lag(cumsum(mean19), default = 0) + mean19) %>% ungroup() p19 <- ggplot(df19, aes(x = Salinity, y = mean19, fill = `Metamorphic response`)) + geom_bar(stat = 'identity') + geom_errorbar(data = error_bars19, aes(ymax = mean19 + se19, ymin = mean19 - se19), position = position_dodge(width = 0.3), width = 0.2) + ggtitle("","") + labs(x = "Salinity (‰)", y = "", fill = NULL) + theme(plot.title = element_text(hjust = 0, size = 24, family = "Times New Roman", face = "bold"), plot.subtitle = element_text(hjust = 0.5, size = 24, family = "Times New Roman"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), axis.text = element_text(size = 30, family = "Times New Roman"), axis.title.y = element_text(size = 30, family = "Times New Roman"), axis.title.x = element_text(size = 30, family = "Times New Roman"), legend.position = "") + # Adjust legend position to the right scale_fill_manual(values = c("#d9d9d9","#737373")) + ylim(0, 120) + scale_y_continuous(limits = c(0, 120),breaks = c(0, 20, 40, 60, 80, 100)) # Set y-axis limits p19 Set_Data_long_7d_7d<- Set_Data_long %>% filter(`Low salinity exposure duration (d)`=="7",`Duration in settlement assay` == "7") Set_Data_long_7d_7d$`Metamorphic response`<-as.factor(Set_Data_long_7d_7d$`Metamorphic response`) df20 <- Set_Data_long_7d_7d %>% group_by(Salinity, `Metamorphic response`) %>% summarise(mean20 = mean(`Value (%)`), se20 = sd(`Value (%)`) / sqrt(n()), # Calculate standard error .groups = 'drop') %>% mutate(Salinity = factor(Salinity, levels = c(25, 30, 34)), `Metamorphic response` = factor(`Metamorphic response`)) df20$`Metamorphic response` <- factor(df20$`Metamorphic response`, levels = c("Metamorphosis", "Juveniles")) # Calculate the position for error bars error_bars20 = df20 %>% arrange(Salinity, desc(`Metamorphic response`)) %>% group_by(Salinity) %>% mutate(mean20 = lag(cumsum(mean20), default = 0) + mean20) %>% ungroup() p20.1 <- ggplot(df20, aes(x = Salinity, y = mean20, fill = `Metamorphic response`)) + geom_bar(stat = 'identity') + geom_errorbar(data = error_bars20, aes(ymax = mean20 + se20, ymin = mean20 - se20), position = position_dodge(width = 0.3), width = 0.2) + ggtitle("","") + labs(x = "", y = "", fill = NULL) + theme(plot.title = element_text(hjust = 0, size = 24, family = "Times New Roman", face = "bold"), plot.subtitle = element_text(hjust = 0.5, size = 24, family = "Times New Roman"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), axis.text = element_text(size = 30, family = "Times New Roman"), axis.title.y = element_text(size = 30, family = "Times New Roman"), axis.title.x = element_text(size = 30, family = "Times New Roman"), legend.position = "") + # Adjust legend position to the right scale_fill_manual(values = c("#d9d9d9","#737373")) + ylim(0, 120) + scale_y_continuous(limits = c(0, 120),breaks = c(0, 20, 40, 60, 80, 100)) # Set y-axis limits p20.1 stackedbps_7d<-plot_grid(p18,p19,p20.1, ncol = 3,nrow = 1) stackedbps_7d<-grid.arrange(arrangeGrob(stackedbps_7d)) ``` ###### Combined plots ```{r} stackedbps<-plot_grid(stackedbps_2d, stackedbps_4d, stackedbps_7d, ncol = 1,nrow = 3) stackedbps ``` #### GLMM ##### Data processing for analyses ```{r} # Transform data columns to factors Set_Data$Sal <- as.factor(Set_Data$`Larval Salinity Group`) Set_Data$Time_Sal <- as.factor(Set_Data$`Pre-settlement assay salinity exposure duration (days)`) Set_Data$Time_SetAssay <- as.factor(Set_Data$`Duration in settlement assay (days)`) Set_Data$cov <- Set_Data$`Handpicked larvae (#normal competent larvae placed in settlement assay)` ID<-as.factor(Set_Data$`Rep ID`) ``` ###### Figure ```{r} # ====== # Fig 6a # ====== y_sett<-(cbind((Set_Data$`Attaching and Resorbing`+Set_Data$Juveniles), Set_Data$`Handpicked larvae (#normal competent larvae placed in settlement assay)`- (Set_Data$`Attaching and Resorbing`+Set_Data$Juveniles))) # Specify the optimization method glmer_control <- glmerControl(optimizer = "bobyqa") #(gold standard for use with repeated measures or convergence issues) #run the rptd measures GLMM #Mod20<-glmer(y_sett ~ Sal*Time_Sal*Time_SetAssay+cov+(1|ID), family = binomial(link = "logit"), data = Set_Data, nAGQ = 25, control = glmer_control) #Mod20<-glmer(y_sett ~ (Sal+Time_Sal+Time_SetAssay)^2+(1|ID), family = binomial(link = "logit"), data = Set_Data, nAGQ = 25, control = glmer_control) Mod20<-glmer(y_sett ~ Sal+Time_Sal+Time_SetAssay+(1|ID), family = binomial(link = "logit"), data = Set_Data, nAGQ = 25, control = glmer_control) Aov20<-Anova(Mod20) Aov20#full interaction term NS emmeans_Mod20<-emmeans(Mod20, ~Time_SetAssay, type = "response") emmeans_Mod20<-emmeans(Mod20, ~Time_Sal, type = "response") emmeans_Mod20<-emmeans(Mod20, ~Sal, type = "response") pairs <- pairs(emmeans_Mod20) pairs emmeans_Mod20<-emmeans(Mod20, ~Time_SetAssay+Sal+Time_Sal, type = "response") emmip(Mod20, ~Time_SetAssay|Sal|Time_Sal, CI=TRUE, type = "response") noperfect_separation<-as.data.frame(emmeans_Mod20) Probs_sett <- ggplot(noperfect_separation, aes(x = Time_SetAssay, y = prob, group = `Sal`, color = `Sal`, fill = `Sal`)) + geom_line(position = position_dodge(width = 0.25)) + geom_point(position = position_dodge(width = 0.25)) + geom_ribbon(aes(ymin = asymp.LCL, ymax = asymp.UCL), alpha = 0.2, position = position_dodge(width = 0.25)) + ggtitle("", "(a)") + ylab("Prob(Metamorphic response)") + xlab("Duration in settlement assay (days)") + facet_wrap(~Time_Sal) + theme(plot.title = element_text(hjust = 0.5, size = 12, family = "Times New Roman", face = "bold"), plot.subtitle = element_text(size = 12, hjust = -0.05, vjust = 0, family = "Times New Roman", face = "bold"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), axis.text = element_text(size = 12, family = "Times New Roman"), strip.background = element_blank(), strip.text = element_text(size = 12, family = "Times New Roman", face = "bold"), axis.title.y = element_text(size = 12, family = "Times New Roman"), axis.title.x = element_text(size = 12, family = "Times New Roman"), legend.position = "right", legend.text = element_text(size = 12, family = "Times New Roman", colour = "black"), legend.title = element_text(size = 12, family = "Times New Roman", colour = "black"), plot.margin = margin(t = 11 * 1.2, r = 5.5, b = 5.5, l = 5.5, unit = "mm")) + # Add top margin scale_color_manual(values = c("#d9d9d9", "#969696", "#252525")) + scale_fill_manual(values = c("#d9d9d9", "#969696", "#252525")) + ylim(0, 1) + labs(color = "Larval duration\nin salinity\ntreatment (days)", fill = "Larval duration\nin salinity\ntreatment (days)") Probs_sett greyscale_colors <- c("#f0f0f0", "#d9d9d9", "#bdbdbd", "#969696", "#737373", "#525252", "#252525") #Assumptions # Fit a model without random effects to check VIF mod_no_re<-glm(y_sett ~ Sal+(Time_Sal+Time_SetAssay)^2, family = binomial(link = "logit"), data = Set_Data) # Calculate VIF for the fixed effects vif_values <- vif(mod_no_re) print(vif_values) r2_mod20 <- r.squaredGLMM(Mod20) print(r2_mod20) # Calculate overdispersion statistic overdispersion_stat <- sum(residuals(Mod20, type = "pearson")^2) / df.residual(Mod20) print(overdispersion_stat) # Residuals vs. fitted values plot(fitted(Mod20), residuals(Mod20, type = "pearson"), xlab = "Fitted values", ylab = "Pearson residuals", main = "Residuals vs Fitted") abline(h = 0, col = "red") ranef_values <- ranef(Mod20, condVar = TRUE) qqnorm(ranef_values$ID[[1]]) qqline(ranef_values$ID[[1]], col = "red") # ====== # Fig 6b # ====== emmeans_Mod20<-emmeans(Mod20, ~Sal, type = "response") emmip(Mod20, ~Sal, CI=TRUE, type = "response") noperfect_separation<-as.data.frame(emmeans_Mod20) Fig_6b <- ggplot(noperfect_separation, aes(x = Sal, y = prob, fill = Sal)) + geom_bar(stat = "identity", position = position_dodge(width = 0.25)) + geom_errorbar(aes(ymin = prob - SE, ymax = prob + SE), width = 0.2, position = position_dodge(0.25)) + ylab("Prob(Metamorphic response)") + xlab("Larval salinity\ntreatment (‰)") + ggtitle("", "(b)") + theme(plot.title = element_text(hjust = 0.5, size = 12, family = "Times New Roman", face = "bold"), plot.subtitle = element_text(size = 40, hjust = -0.05, vjust = 5, family = "Times New Roman", face = "bold"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), axis.text = element_text(size = 40, family = "Times New Roman"), strip.background = element_blank(), strip.text = element_text(size = 40, family = "Times New Roman", face = "bold"), axis.title.y = element_text(size = 40, family = "Times New Roman"), axis.title.x = element_text(size = 40, family = "Times New Roman"), legend.position = "", legend.text = element_text(size = 40, family = "Times New Roman"), legend.title = element_text(size = 40, family = "Times New Roman")) + ylim(0, 1) + scale_fill_manual(values = c("grey", "grey", "grey")) + scale_colour_manual(values = c("grey", "grey", "grey")) + annotate("text", x = 1, y = 0.33, label = "a", size = 15, family = "Times New Roman") + annotate("text", x = 2, y = 0.55, label = "b", size = 15, family = "Times New Roman") + annotate("text", x = 3, y = 0.48, label = "ab", size = 15, family = "Times New Roman") Fig_6b # ====== # Fig 6c # ====== emmeans_Mod20<-emmeans(Mod20, ~Time_Sal, type = "response") emmip(Mod20, ~Time_Sal, CI=TRUE, type = "response") noperfect_separation<-as.data.frame(emmeans_Mod20) Fig_6c <- ggplot(noperfect_separation, aes(x = Time_Sal, y = prob, fill = Time_Sal)) + geom_bar(stat = "identity", position = position_dodge(width = 0.25)) + geom_errorbar(aes(ymin = prob - SE, ymax = prob + SE), width = 0.2, position = position_dodge(0.25)) + ylab("Prob(Metamorphic response)") + xlab("Larval duration in\nsalinity treatment (days)") + ggtitle("", "(c)") + theme(plot.title = element_text(hjust = 0.5, size = 12, family = "Times New Roman", face = "bold"), plot.subtitle = element_text(size = 40, hjust = -0.05, vjust = 5, family = "Times New Roman", face = "bold"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), axis.text = element_text(size = 40, family = "Times New Roman"), strip.background = element_blank(), strip.text = element_text(size = 40, family = "Times New Roman", face = "bold"), axis.title.y = element_text(size = 40, family = "Times New Roman"), axis.title.x = element_text(size = 40, family = "Times New Roman"), legend.position = "", legend.text = element_text(size = 40, family = "Times New Roman"), legend.title = element_text(size = 40, family = "Times New Roman")) + ylim(0, 1) + scale_fill_manual(values = c("grey", "grey", "grey")) + scale_colour_manual(values = c("grey", "grey", "grey")) + annotate("text", x = 1, y = 0.30, label = "a", size = 15, family = "Times New Roman") + annotate("text", x = 2, y = 0.53, label = "b", size = 15, family = "Times New Roman") + annotate("text", x = 3, y = 0.55, label = "b", size = 15, family = "Times New Roman") Fig_6c # ====== # Fig 6d # ====== emmeans_Mod20<-emmeans(Mod20, ~Time_SetAssay, type = "response") emmip(Mod20, ~Time_SetAssay, CI=TRUE, type = "response") noperfect_separation<-as.data.frame(emmeans_Mod20) Fig_6d <- ggplot(noperfect_separation, aes(x = Time_SetAssay, y = prob, fill = Time_SetAssay)) + geom_bar(stat = "identity", position = position_dodge(width = 0.25)) + geom_errorbar(aes(ymin = prob - SE, ymax = prob + SE), width = 0.2, position = position_dodge(0.25)) + ylab("Prob(Metamorphic response)") + xlab("Duration in\nsettlement assays (days)") + ggtitle("", "(d)") + theme(plot.title = element_text(hjust = 0.5, size = 12, family = "Times New Roman", face = "bold"), plot.subtitle = element_text(size = 40, hjust = -0.05, vjust = 5, family = "Times New Roman", face = "bold"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), axis.text = element_text(size = 40, family = "Times New Roman"), strip.background = element_blank(), strip.text = element_text(size = 40, family = "Times New Roman", face = "bold"), axis.title.y = element_text(size = 40, family = "Times New Roman"), axis.title.x = element_text(size = 40, family = "Times New Roman"), legend.position = "", legend.text = element_text(size = 40, family = "Times New Roman"), legend.title = element_text(size = 40, family = "Times New Roman")) + ylim(0, 1) + scale_fill_manual(values = c("grey", "grey", "grey")) + scale_colour_manual(values = c("grey", "grey", "grey")) + annotate("text", x = 1, y = 0.32, label = "a", size = 15, family = "Times New Roman") + annotate("text", x = 2, y = 0.46, label = "b", size = 15, family = "Times New Roman") + annotate("text", x = 3, y = 0.53, label = "b", size = 15, family = "Times New Roman") Fig_6d ``` ### Juvenile success #### Load data ```{r} # Load the data Juv_Data <- read_excel('[insert filepath]/Clements & Byrne 2026.xlsx', sheet = "Juvenile abundance", range = "A2:F130", col_names = c("Rep ID", "Pre-settlement assay salinity exposure duration (days)", "Duration in settlement assay (days)", "Larval Salinity Group", "Handpicked larvae (#normal competent larvae placed in settlement assay)", "Juveniles")) # Transform data columns to factors Juv_Data$Salinity <- as.factor(Juv_Data$`Larval Salinity Group`) Juv_Data$Time_exposed <- as.factor(Juv_Data$`Pre-settlement assay salinity exposure duration (days)`) Juv_Data$Time_in_assay <- as.factor(Juv_Data$`Duration in settlement assay (days)`) Juv_Data$cov <- Juv_Data$`Handpicked larvae (#normal competent larvae placed in settlement assay)` ``` #### GLMM ```{r} # ======== # Table 4b # ======== y_juv<-(cbind((Juv_Data$Juveniles), Juv_Data$`Handpicked larvae (#normal competent larvae placed in settlement assay)`- (Juv_Data$Juveniles))) Sal<-as.factor(Juv_Data$`Larval Salinity Group`) Time_Sal<-as.factor(Juv_Data$`Pre-settlement assay salinity exposure duration (days)`) Time_sett<-as.factor(Juv_Data$`Duration in settlement assay (days)`) cov<-Juv_Data$`Handpicked larvae (#normal competent larvae placed in settlement assay)` ID<-as.factor(Juv_Data$`Rep ID`) Mod9<-glmer(y_juv ~ Sal*Time_Sal*Time_sett+cov+(1|ID), family = binomial(link = "logit"), data = Set_Data, nAGQ = 25, control = glmer_control) Mod9<-glmer(y_juv ~ (Sal+Time_Sal+Time_sett)^2+(1|ID), family = binomial(link = "logit"), data = Set_Data, nAGQ = 25, control = glmer_control) Mod9<-glmer(y_juv ~ (Sal+Time_Sal)^2+(Time_Sal+Time_sett)^2+(1|ID), family = binomial(link = "logit"), data = Juv_Data, nAGQ = 25, control = glmer_control) Aov9<-Anova(Mod9) Aov9 emmeans_Mod9<-emmeans(Mod9, ~ (Sal+Time_Sal)^2+(Time_Sal+Time_sett)^2, type = "response") emmeans_Mod9<-emmeans(Mod9, ~Sal, type = "response") emmeans_Mod9<-emmeans(Mod9, ~Time_Sal, type = "response") emmeans_Mod9<-emmeans(Mod9, ~Time_sett, type = "response") pairs <- pairs(emmeans_Mod9) pairs #emmip(Mod9, ~ Sal|Time_Sal|Time_sett, CI=TRUE, type = "response") #Perfect seperation detected # ========= # Figure 7a # ========= emmeans_Mod9<-emmeans(Mod9, ~ (Sal+Time_Sal)^2+(Time_Sal+Time_sett)^2, type = "response") noperfect_separation_juv<-as.data.frame(emmeans_Mod9) Figure_7a <- ggplot(noperfect_separation_juv, aes(x = Time_sett, y = prob, group = `Sal`, color = `Sal`)) + geom_line(position = position_dodge(width = 0.25)) + geom_point(position = position_dodge(width = 0.25)) + geom_ribbon(aes(ymin = asymp.LCL, ymax = asymp.UCL, fill = `Sal`), alpha = 0.2, position = position_dodge(width = 0.25)) + ylab("Prob(Juvenile)") + xlab("Duration in settlement assay (days)") + ggtitle("", "(a)") + facet_wrap(~Time_Sal)+ theme(plot.title = element_text(hjust = 0.5, size = 12, family = "Times New Roman", face = "bold"), plot.subtitle = element_text(size = 12, hjust = -0.05, vjust = -10, family = "Times New Roman", face = "bold"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), axis.text = element_text(size = 12, family = "Times New Roman"), strip.background = element_blank(), strip.text = element_text(size = 12, family = "Times New Roman"), axis.title.y = element_text(size = 12, family = "Times New Roman"), axis.title.x = element_text(size = 12, family = "Times New Roman"), legend.position = "right", legend.text = element_text(size = 12, family = "Times New Roman"), legend.title = element_text(size = 12, family = "Times New Roman")) + scale_fill_manual(values = c("#d9d9d9", "#969696", "#252525"), name = "Larval salinity\ntreatment (‰)") + scale_colour_manual(values = c("#d9d9d9", "#969696", "#252525"), name = "Larval salinity\ntreatment (‰)") + ylim(0, 1) Figure_7a # ========= # Figure 7b # ========= emmeans_Mod9<-emmeans(Mod9, ~Sal, type = "response") noperfect_separation_juv<-as.data.frame(emmeans_Mod9) Fig_7b <- ggplot(noperfect_separation_juv, aes(x = Sal, y = prob, fill = Sal)) + geom_bar(stat = "identity", position = position_dodge(width = 0.25)) + geom_errorbar(aes(ymin = prob - SE, ymax = prob + SE), width = 0.2, position = position_dodge(0.25)) + ylab("Prob(Juvenile)") + xlab("Larval salinity\ntreatment (‰)") + ggtitle("", "(b)") + theme(plot.title = element_text(hjust = 0.5, size = 12, family = "Times New Roman", face = "bold"), plot.subtitle = element_text(size = 40, hjust = -0.05, vjust = 5, family = "Times New Roman", face = "bold"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), axis.text = element_text(size = 40, family = "Times New Roman"), strip.background = element_blank(), strip.text = element_text(size = 40, family = "Times New Roman", face = "bold"), axis.title.y = element_text(size = 40, family = "Times New Roman"), axis.title.x = element_text(size = 40, family = "Times New Roman"), legend.position = "", legend.text = element_text(size = 40, family = "Times New Roman"), legend.title = element_text(size = 40, family = "Times New Roman")) + ylim(0, 1) + scale_fill_manual(values = c("grey", "grey", "grey")) + scale_colour_manual(values = c("grey", "grey", "grey")) + annotate("text", x = 1, y = 0.18, label = "a", size = 15, family = "Times New Roman") + annotate("text", x = 2, y = 0.29, label = "b", size = 15, family = "Times New Roman") + annotate("text", x = 3, y = 0.25, label = "b", size = 15, family = "Times New Roman") Fig_7b # ========= # Figure 7c # ========= emmeans_Mod9<-emmeans(Mod9, ~Time_Sal, type = "response") noperfect_separation_juv<-as.data.frame(emmeans_Mod9) Fig_7c <- ggplot(noperfect_separation_juv, aes(x = Time_Sal, y = prob, fill = Time_Sal)) + geom_bar(stat = "identity", position = position_dodge(width = 0.25)) + geom_errorbar(aes(ymin = prob - SE, ymax = prob + SE), width = 0.2, position = position_dodge(0.25)) + ylab("Prob(Juvenile)") + xlab("Larval duration\nin salinity treatment (days)") + ggtitle("", "(c)") + theme(plot.title = element_text(hjust = 0.5, size = 12, family = "Times New Roman", face = "bold"), plot.subtitle = element_text(size = 40, hjust = -0.05, vjust = 5, family = "Times New Roman", face = "bold"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), axis.text = element_text(size = 40, family = "Times New Roman"), strip.background = element_blank(), strip.text = element_text(size = 40, family = "Times New Roman", face = "bold"), axis.title.y = element_text(size = 40, family = "Times New Roman"), axis.title.x = element_text(size = 40, family = "Times New Roman"), legend.position = "", legend.text = element_text(size = 40, family = "Times New Roman"), legend.title = element_text(size = 40, family = "Times New Roman")) + ylim(0, 1) + scale_fill_manual(values = c("grey", "grey", "grey")) + scale_colour_manual(values = c("grey", "grey", "grey")) + annotate("text", x = 1, y = 0.15, label = "a", size = 15, family = "Times New Roman") + annotate("text", x = 2, y = 0.20, label = "b", size = 15, family = "Times New Roman") + annotate("text", x = 3, y = 0.30, label = "b", size = 15, family = "Times New Roman") Fig_7c # ========= # Figure 7d # ========= emmeans_Mod9<-emmeans(Mod9, ~Time_sett, type = "response") noperfect_separation_juv<-as.data.frame(emmeans_Mod9) Fig_7d <- ggplot(noperfect_separation_juv, aes(x = Time_sett, y = prob, fill = Time_sett)) + geom_bar(stat = "identity", position = position_dodge(width = 0.25)) + geom_errorbar(aes(ymin = prob - SE, ymax = prob + SE), width = 0.2, position = position_dodge(0.25)) + ylab("Prob(Juvenile)") + xlab("Duration\nin settlement assay (days)") + ggtitle("", "(d)") + theme(plot.title = element_text(hjust = 0.5, size = 12, family = "Times New Roman", face = "bold"), plot.subtitle = element_text(size = 40, hjust = -0.05, vjust = 5, family = "Times New Roman", face = "bold"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), axis.text = element_text(size = 40, family = "Times New Roman"), strip.background = element_blank(), strip.text = element_text(size = 40, family = "Times New Roman", face = "bold"), axis.title.y = element_text(size = 40, family = "Times New Roman"), axis.title.x = element_text(size = 40, family = "Times New Roman"), legend.position = "", legend.text = element_text(size = 40, family = "Times New Roman"), legend.title = element_text(size = 40, family = "Times New Roman")) + ylim(0, 1) + scale_fill_manual(values = c("grey", "grey", "grey")) + scale_colour_manual(values = c("grey", "grey", "grey")) + annotate("text", x = 1, y = 0.10, label = "a", size = 15, family = "Times New Roman") + annotate("text", x = 2, y = 0.28, label = "b", size = 15, family = "Times New Roman") + annotate("text", x = 3, y = 0.35, label = "c", size = 15, family = "Times New Roman") Fig_7d ``` # ==================== # Juvenile production # ==================== ### 1d old juvies ### Mean duration to reach the juvenile stage ```{r} Juv_Data_COE1 <- read_excel('[insert filepath]/Clements & Byrne 2026.xlsx', sheet = "Juvenile production & size", range = "A2:G169", col_names = c("Pre-settlement assay salinity exposure duration (days)", "Elapsed time to become a juvenile (duration in settlement assay (days))", "Larval Salinity Group", "Rep unique code", "Handpicked", "Days since settlement", "Area (10^5 um^2)")) Juv_Data <- Juv_Data_COE1 %>% filter(`Days since settlement` == "1 day old" & `Larval Salinity Group` %in% c("30", "34")) Juv_Data$Salinity <- as.factor(Juv_Data$`Larval Salinity Group`) Juv_Data$Time_exposed <- as.factor(Juv_Data$`Pre-settlement assay salinity exposure duration (days)`) Juv_Data$elapsedtime <- Juv_Data$`Elapsed time to become a juvenile (duration in settlement assay (days))` # ===================== # Supplementary Table 6 # ===================== # Calculate mean elapsed time to become a juvenile by salinity group mean_elapsed_time <- Juv_Data %>% group_by(Salinity, `Pre-settlement assay salinity exposure duration (days)`) %>% summarize( mean_elapsed_time = mean(elapsedtime, na.rm = TRUE), std_error = sd(elapsedtime, na.rm = TRUE) / sqrt(n()), sample_size = n() ) # ===================================== # ANCOVA: duration to become a juvenile # ===================================== anova_result1 <- aov(elapsedtime ~ (Salinity * Time_exposed) + `Handpicked`, data = Juv_Data) summary(anova_result1) par(mfrow = c(2, 2)) plot(anova_result1) Juv_Data$resid_anova1 <- residuals(anova_result1) leveneTest(resid_anova1 ~ Salinity * Time_exposed, data = Juv_Data) # Perform pairwise comparisons emmeans_result <- emmeans(anova_result1, ~ Salinity * Time_exposed) # Pairwise salinity comparisons within each exposure time pairs(emmeans(anova_result1, ~ Salinity | Time_exposed)) # Pairwise exposure-time comparisons within each salinity pairs(emmeans(anova_result1, ~ Time_exposed | Salinity)) ``` ## Plot ```{r} # Ensure the factor levels are set correctly in the data mean_elapsed_time$`Pre-settlement assay salinity exposure duration (days)` <- factor(mean_elapsed_time$`Pre-settlement assay salinity exposure duration (days)`, levels = c("2d", "4d", "7d")) Duration_to_juvie <- ggplot(mean_elapsed_time, aes(x = Salinity, y = mean_elapsed_time, fill = factor(`Pre-settlement assay salinity exposure duration (days)`))) + geom_bar(stat = "identity", position = position_dodge(), width = 0.7) + geom_errorbar(aes(ymin = mean_elapsed_time - std_error, ymax = mean_elapsed_time + std_error), position = position_dodge(width = 0.7), width = 0.25) + labs(title = "", x = "Larval salinity treatment (PSU)", y = "Time until juvenile (days)", fill = "Larval duration in\nsalinity treatment") + scale_fill_manual(values = c("2d" = "#d9d9d9", "4d" = "#969696", "7d" = "#525252"), labels = c("2 days", "4 days", "7 days")) + scale_y_continuous(limits = c(0, 6), breaks = seq(0, 7, by = 1)) + theme(plot.title = element_text(hjust = 0, size = 12, family = "Times New Roman", face = "bold"), plot.subtitle = element_text(hjust = 0.5, size = 12, family = "Times New Roman"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), axis.text = element_text(size = 12, family = "Times New Roman"), axis.title.y = element_text(size = 12, family = "Times New Roman"), axis.title.x = element_text(size = 12, family = "Times New Roman"), legend.title = element_text(size = 12, family = "Times New Roman"), legend.text = element_text(size = 12, family = "Times New Roman"), legend.position = "right") Duration_to_juvie summary_data <- Juv_Data %>% group_by(Salinity, Time_exposed) %>% summarise(mean_elapsed_time = mean(`Elapsed time to become a juvenile (duration in settlement assay (days))`), se_elapsed_time = sd(`Elapsed time to become a juvenile (duration in settlement assay (days))`)/sqrt(n())) ``` # ============= # Juvenile size # ============= ### 1 day old juveniles ```{r} # ======== # Table 6a # ======== oneday_old <- Juv_Data_COE1 %>% filter(`Pre-settlement assay salinity exposure duration (days)` %in% c("2d", "4d"), `Days since settlement` == "1 day old", `Larval Salinity Group` %in% c("30", "34")) # Transform relevant columns to factors oneday_old$Salinity <- as.factor(oneday_old$`Larval Salinity Group`) oneday_old$Time_exposed <- as.factor(oneday_old$`Pre-settlement assay salinity exposure duration (days)`) oneday_old$Age<-as.factor(oneday_old$`Days since settlement`) oneday_old$Rep<-as.factor(oneday_old$`Rep unique code`) # Calculate mean mean_Area_day1 <- oneday_old %>% group_by(Salinity, Time_exposed, Age) %>% summarize(mean_area1day = mean(`Area (10^5 um^2)`, na.rm = TRUE), std_error = sd(`Area (10^5 um^2)`, na.rm = TRUE) / sqrt(n()), sample_size = n()) mean_Area_day1 # Fit the model (if not already done) lme_result <- lmer(`Area (10^5 um^2)` ~ Salinity * Time_exposed + (1|Rep) , data = oneday_old) lme_result <- lmer(`Area (10^5 um^2)` ~ Salinity + Time_exposed + (1|Rep), data = oneday_old) # Perform an ANOVA-like table of fixed effects Anova(lme_result) # Residuals vs Fitted plot plot(fitted(lme_result), residuals(lme_result)) abline(h = 0, col = "red") shapiro.test(residuals(lme_result)) Area_day1 <- ggplot(mean_Area_day1, aes(x = Salinity, y = mean_area1day, fill = Time_exposed)) + geom_bar(stat = "identity", position = position_dodge(), width = 0.7) + geom_errorbar(aes(ymin = mean_area1day - std_error, ymax = mean_area1day + std_error), position = position_dodge(width = 0.7), width = 0.25) + labs(title = "(a) 1 day old", x = "Larval salinity treatment (PSU)", y = expression("Aboral body area (10"^5 ~ mu*m^2 ~ ")"), fill = "Larval salinity treatment \nduration (days)") + scale_fill_manual(values = c("#d9d9d9","#969696")) + theme(plot.title = element_text(hjust = 0, size = 12, family = "Times New Roman", face = "bold"), plot.subtitle = element_text(hjust = 0.5, size = 12, family = "Times New Roman"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), axis.text = element_text(size = 12, family = "Times New Roman"), axis.title.y = element_text(size = 12, family = "Times New Roman"), axis.title.x = element_text(size = 12, family = "Times New Roman"), strip.background = element_blank(), strip.text = element_blank(), legend.text = element_text(size = 12, family = "Times New Roman"), legend.title = element_text(size = 12, family = "Times New Roman"), legend.position = "")+ ylim(0, 4) # Print the plot print(Area_day1) ``` ### 5 day old juveniles ```{r} # Load the data from the specified Excel file Juv_Data_COE5 <- read_excel('[insert filepath]/Clements & Byrne 2026.xlsx', sheet = "Juvenile production & size", range = "A58:G111", col_names = c("Pre-settlement assay salinity exposure duration (days)", "Elapsed time to become a juvenile (duration in settlement assay (days))", "Larval Salinity Group", "Rep unique code", "Handpicked", "Days since settlement", "Area (10^5 um^2)")) fivedays_old <- Juv_Data_COE5 %>% filter(`Pre-settlement assay salinity exposure duration (days)` %in% c("2d", "4d"), `Larval Salinity Group` %in% c("30", "34")) # Transform relevant columns to factors fivedays_old$Salinity <- as.factor(fivedays_old$`Larval Salinity Group`) fivedays_old$Time_exposed <- as.factor(fivedays_old$`Pre-settlement assay salinity exposure duration (days)`) fivedays_old$Age<-as.factor(fivedays_old$`Days since settlement`) fivedays_old$Rep<-as.factor(fivedays_old$`Rep unique code`) # Calculate mean mean_Area_day5 <- fivedays_old %>% group_by(Salinity, Time_exposed, Age) %>% summarize(mean_area5day = mean(`Area (10^5 um^2)`, na.rm = TRUE), std_error = sd(`Area (10^5 um^2)`, na.rm = TRUE) / sqrt(n()), sample_size = n()) mean_Area_day5 # Fit the model (if not already done) lme_result <- lmer(`Area (10^5 um^2)` ~ Salinity * Time_exposed + (1|Rep), data = fivedays_old) lme_result <- lmer(`Area (10^5 um^2)` ~ Salinity + Time_exposed + (1|Rep), data = fivedays_old) # Perform an ANOVA-like table of fixed effects Anova(lme_result) anova_result4<-Anova(lme_result) emmeans<-emmeans(lme_result, ~Salinity, adjust = "tukey") pairs <- pairs(emmeans) pairs # Residuals vs Fitted plot plot(fitted(lme_result), residuals(lme_result)) abline(h = 0, col = "red") qqnorm(residuals(lme_result)) qqline(residuals(lme_result), col = "red") #Normality shapiro.test(residuals(lme_result)) Area_day5 <- ggplot(mean_Area_day5, aes(x = Salinity, y = mean_area5day, fill = Time_exposed)) + geom_bar(stat = "identity", position = position_dodge(), width = 0.7) + geom_errorbar(aes(ymin = mean_area5day - std_error, ymax = mean_area5day + std_error), position = position_dodge(width = 0.7), width = 0.25) + labs(title = "(b) 5 days old", x = "Larval salinity treatment (‰)", y = expression(""), fill = "Larval salinity\ntreatment duration") + scale_fill_manual(values = c("#d9d9d9","#969696")) + theme(plot.title = element_text(hjust = 0, size = 12, family = "Times New Roman", face = "bold"), plot.subtitle = element_text(hjust = 0.5, size = 12, family = "Times New Roman"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"), axis.text = element_text(size = 12, family = "Times New Roman"), axis.title.y = element_text(size = 12, family = "Times New Roman"), axis.title.x = element_text(size = 12, family = "Times New Roman"), strip.background = element_blank(), strip.text = element_blank(), legend.text = element_text(size = 12, family = "Times New Roman"), legend.title = element_text(size = 12, family = "Times New Roman"), legend.position = "") # Print the plot print(Area_day5) ``` #### Combined plots ```{r} Bodysize <- plot_grid(Area_day1, Area_day5, ncol = 2, rel_widths = c(1, 1)) Bodysize ```