---
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
```
