---
title: Teratogenic effects of larval low salinity experience in development of the
  juvenile body in Acanthaster sp.
author: "Clements et al., 2025"
date: "`r Sys.Date()`"
output: html_document

#Description
##The code in the chunks below are provided to support open and reproducable science. The code was authored in full by Matthew Clements, the lead author of the manuscript titled - "Teratogenic effects of larval low salinity experience in development of the juvenile body in Acanthaster sp.". The code in this file supports the data analysed and visualised within this manuscript. The data will be available in the USYD (Sydney eScholarship Repository). Please locate the published manuscript for the DOI to access the data or contact matthew.clements@sydney.edu.au.   

---
#Load required packages
```{r}
library(readxl)
library(dplyr)
library(ggplot2)
library(cowplot)
library(lme4)
library(car)
library(emmeans)
library(DHARMa)
library(lattice)
library(arm)
```

#Influence of the larval salinity experience on juvenile growth
##1 day old juveniles
##Load data 
```{r}
Data_1day_old <- read_excel("/Users/mcle5454/Desktop/For submission/Brach project 2/Post revision/Teratogenic effects of larval low salinity experience in development of the juvenile body in Acanthaster sp..xlsx", 
                           sheet = "Juvenile growth", 
                           range = "A2:H104", 
                           col_names = c("Ref #",	"Juvie Id",	"Larval salinity treatment",	"Larval salinity exposure duration",	"Replicate#",	"Juvenile#",	"Juvenile age (days)",	"Body diameter (arm tip-arm tip, um)"))

# Transform relevant columns to factors
Data_1day_old$`Larval salinity treatment` <- as.factor(Data_1day_old$`Larval salinity treatment`)
Data_1day_old$`Larval salinity exposure duration` <- as.factor(Data_1day_old$`Larval salinity exposure duration`)
```

###Analyses
####LMM_1day_old (body diameter)
```{r}
Data_1day_old$sal<-Data_1day_old$`Larval salinity treatment`
Data_1day_old$time<-Data_1day_old$`Larval salinity exposure duration`

lmm_model1 <- lmer(`Body diameter (arm tip-arm tip, um)` ~ sal * time + (1 | `Replicate#`), data = Data_1day_old) #full model (i.e. not the selected model: y~salinity*time+(1|ID)). 
Aov1<-Anova(lmm_model1)
Aov1 #The interaction term was dropped for the best model fit (lmm_model2)

lmm_model2 <- lmer(`Body diameter (arm tip-arm tip, um)` ~ sal + time + (1 | `Replicate#`), data = Data_1day_old) #best model fit (i.e: y~salinity+time+(1|ID))
Aov2<-Anova(lmm_model2)
Aov2

emmip(lmm_model2, ~sal|time, CI=TRUE, type = "response") 

emmeans_Mod2<-emmeans(lmm_model2, ~sal, type = "response") 
emmeans_Mod2<-emmeans(lmm_model2, ~time, type = "response") 

pairs <- pairs(emmeans_Mod2)
pairs

noperfect_separation_juv_1<-as.data.frame(emmeans_Mod1)
```
#####Assumptions
```{r}
#The following code tests the assumptions of the above analysis titled "LMM_1day_old (body diameter)":
##Homoscedasticity (constant variance of residuals)
##Normality of residuals
##Multicollinearity among fixed effects
##Model fit and residual diagnostics

# 1. Residual vs. Fitted plot (homoscedasticity check)
residuals <- resid(lmm_model2)
fitted_values <- fitted(lmm_model2)

ggplot(data = data.frame(Fitted = fitted_values, Residuals = residuals), aes(x = Fitted, y = Residuals)) +
  geom_point() +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  labs(title = "Residuals vs. Fitted", x = "Fitted Values", y = "Residuals") +
  theme_minimal()

# 2. QQ Plot (normality of residuals)
qqnorm(residuals)
qqline(residuals, col = "red")

# 3. Multicollinearity check (VIF for fixed effects)
# Create a linear model for VIF calculation (fixed effects only)
vif_model <- vif(lm(`Body diameter (arm tip-arm tip, um)` ~ sal + time, data = Data_1day_old))
print(vif_model)

# 4. Homoscedasticity check (Scale-Location Plot)
ggplot(data = data.frame(Fitted = fitted_values, Residuals = abs(residuals)), 
       aes(x = Fitted, y = Residuals)) +
  geom_point() +
  geom_smooth(method = "loess", color = "red", se = FALSE) +
  labs(title = "Scale-Location Plot", x = "Fitted Values", y = "Absolute Residuals") +
  theme_minimal()

# 5. Residual diagnostics with DHARMa
# Simulate residuals
simulation_output <- simulateResiduals(fittedModel = lmm_model1, n = 1000)

# Diagnostic plots
plot(simulation_output)

# Overdispersion test
testDispersion(simulation_output)

# Uniformity test (Kolmogorov-Smirnov)
testUniformity(simulation_output)
```
###Mean arm diameter_1 day old juveniles 
```{r}
#All 1 day old/2 day exposure data
Data_1day_old_2day_salinity_exposure <-
  Data_1day_old %>% 
  filter(`Larval salinity exposure duration` == "2 days")

Mean_bodydiameter_1day_old_2day_salinity_exposure <- Data_1day_old_2day_salinity_exposure %>% 
  group_by(`Larval salinity treatment`) %>% 
  summarize(
    mean_armdiam1 = mean(`Body diameter (arm tip-arm tip, um)`, na.rm = TRUE),
    std_error = sd(`Body diameter (arm tip-arm tip, um)`, na.rm = TRUE) / sqrt(n()),
    sample_size = n()
  )

print(Mean_bodydiameter_1day_old_2day_salinity_exposure)

#All 1 day old/4 day exposure data
Data_1day_old_4day_salinity_exposure <-
  Data_1day_old %>% 
  filter(`Larval salinity exposure duration` == "4 days")

mean_bodydiameter_1day_old_4day_salinity_exposure <- Data_1day_old_4day_salinity_exposure %>% 
  group_by(`Larval salinity treatment`) %>% 
  summarize(
    mean_armdiam1 = mean(`Body diameter (arm tip-arm tip, um)`, na.rm = TRUE),
    std_error = sd(`Body diameter (arm tip-arm tip, um)`, na.rm = TRUE) / sqrt(n()),
    sample_size = n()
  )

print(mean_bodydiameter_1day_old_4day_salinity_exposure)
```

#####Mean body diameter in each salinity, regardless of exposure time
```{r}
mean_bodydiameter_1day_old_sal <- Data_1day_old %>% 
  group_by(`Larval salinity treatment`) %>% 
  summarize(
    mean_armdiam1 = mean(`Body diameter (arm tip-arm tip, um)`, na.rm = TRUE),
    std_error = sd(`Body diameter (arm tip-arm tip, um)`, na.rm = TRUE) / sqrt(n()),
    sample_size = n()
  )

print(mean_bodydiameter_1day_old_sal)

#34-30ppt
((550.6087-481.1579)/481.1579)*100

#34-22ppt
((550.6087-440.2857)/440.2857)*100
```

#####Mean body diameter by exposure time, regardless of salinity
```{r}
mean_bodydiameter_1day_old_2day_salinity_exposure_time <- Data_1day_old_2day_salinity_exposure %>% 
  group_by(`Larval salinity exposure duration`) %>% 
  summarize(
    mean_armdiam1 = mean(`Body diameter (arm tip-arm tip, um)`, na.rm = TRUE),
    std_error = sd(`Body diameter (arm tip-arm tip, um)`, na.rm = TRUE) / sqrt(n()),
    sample_size = n()
  )

print(mean_bodydiameter_1day_old_2day_salinity_exposure_time)

mean_bodydiameter_1day_old_4day_salinity_exposure_time <- Data_1day_old_4day_salinity_exposure %>% 
  group_by(`Larval salinity exposure duration`) %>% 
  summarize(
    mean_armdiam2 = mean(`Body diameter (arm tip-arm tip, um)`, na.rm = TRUE),
    std_error = sd(`Body diameter (arm tip-arm tip, um)`, na.rm = TRUE) / sqrt(n()),
    sample_size = n()
  )

print(mean_bodydiameter_1day_old_4day_salinity_exposure_time)

#2d-4d
((509.8-463.9)/463.9)*100
509.8-463.9

```

##5 day old juveniles
##Load data 
```{r}
Data_5day_old <- read_excel("/Users/mcle5454/Desktop/For submission/Brach project 2/Post revision/Teratogenic effects of larval low salinity experience in development of the juvenile body in Acanthaster sp..xlsx", 
                           sheet = "Juvenile growth", 
                           range = "A105:H207", 
                           col_names = c("Ref #",	"Juvie Id",	"Larval salinity treatment",	"Larval salinity exposure duration",	"Replicate#",	"Juvenile#",	"Juvenile age (days)",	"Body diameter (arm tip-arm tip, um)"))

# Transform relevant columns to factors
Data_5day_old$`Larval salinity treatment` <- as.factor(Data_5day_old$`Larval salinity treatment`)
Data_5day_old$`Larval salinity exposure duration` <- as.factor(Data_5day_old$`Larval salinity exposure duration`)

# Transform relevant columns to factors
Data_5day_old$`Larval salinity treatment` <- as.factor(Data_5day_old$`Larval salinity treatment`)
Data_5day_old$`Larval salinity exposure duration` <- as.factor(Data_5day_old$`Larval salinity exposure duration`)
Data_5day_old$`Body diameter (arm tip-arm tip, um)` <- as.numeric(Data_5day_old$`Body diameter (arm tip-arm tip, um)`)

# Remove rows with NA in the specified column
Data_5day_old <- Data_5day_old %>%
  filter(!is.na(`Body diameter (arm tip-arm tip, um)`))
```

###LMM_5day_old (body diameter)
```{r}
Data_5day_old <- Data_5day_old %>%
  filter(!is.na(suppressWarnings(as.numeric(`Body diameter (arm tip-arm tip, um)`)))) %>%
  mutate(`Body diameter (arm tip-arm tip, um)` = as.numeric(`Body diameter (arm tip-arm tip, um)`))

Data_5day_old$sal<-Data_5day_old$`Larval salinity treatment`
Data_5day_old$time<-Data_5day_old$`Larval salinity exposure duration`

lmm_model3 <- lmer(`Body diameter (arm tip-arm tip, um)` ~ sal * time + (1 |`Replicate#`), data = Data_5day_old) #full model (i.e. not the selected model: y~salinity*time+(1|ID)).
Aov3<-Anova(lmm_model3)
Aov3 #The interaction term was dropped for the best model fit (lmm_model4)

lmm_model4 <- lmer(`Body diameter (arm tip-arm tip, um)` ~ sal + time + (1 | `Replicate#`), data = Data_5day_old) #best model fit (i.e: y~salinity+time+(1|ID))
Aov4<-Anova(lmm_model4)
Aov4

emmeans_Mod4<-emmeans(lmm_model4, ~sal+time, type = "response") 
emmip(lmm_model4, ~sal|time, CI=TRUE, type = "response") 

emmeans_Mod4.1<-emmeans(lmm_model4, ~sal, type = "response") 
pairs <- pairs(emmeans_Mod4.1)
pairs

emmeans_Mod4.2<-emmeans(lmm_model4, ~time, type = "response")
pairs <- pairs(emmeans_Mod4.2)
pairs

noperfect_separation_juv_5<-as.data.frame(emmeans_Mod4.1)
```

####Assumptions
```{r}
#The following code tests the assumptions of the above analysis titled "LMM_1day_old (body diameter)":
##Homoscedasticity (constant variance of residuals)
##Normality of residuals
##Multicollinearity among fixed effects
##Model fit and residual diagnostics

# Load required libraries
# 1. Residual vs. Fitted plot (homoscedasticity check)
residuals <- resid(lmm_model4)
fitted_values <- fitted(lmm_model4)

ggplot(data = data.frame(Fitted = fitted_values, Residuals = residuals), aes(x = Fitted, y = Residuals)) +
  geom_point() +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  labs(title = "Residuals vs. Fitted", x = "Fitted Values", y = "Residuals") +
  theme_minimal()

# 2. QQ Plot (normality of residuals)
qqnorm(residuals)
qqline(residuals, col = "red")

# 3. Multicollinearity check (VIF for fixed effects)
# Create a linear model for VIF calculation (fixed effects only)
vif_model <- vif(lm(`Body diameter (arm tip-arm tip, um)` ~ sal + time, data = Data_5day_old))
print(vif_model)

# 4. Homoscedasticity check (Scale-Location Plot)
ggplot(data = data.frame(Fitted = fitted_values, Residuals = abs(residuals)), 
       aes(x = Fitted, y = Residuals)) +
  geom_point() +
  geom_smooth(method = "loess", color = "red", se = FALSE) +
  labs(title = "Scale-Location Plot", x = "Fitted Values", y = "Absolute Residuals") +
  theme_minimal()

# 5. Residual diagnostics with DHARMa
# Simulate residuals
simulation_output <- simulateResiduals(fittedModel = lmm_model1, n = 1000)

# Diagnostic plots
plot(simulation_output)

# Overdispersion test
testDispersion(simulation_output)

# Uniformity test (Kolmogorov-Smirnov)
testUniformity(simulation_output)
```

###Mean arm diameter_5 day old juveniles 
```{r}
#All 1 day old/2 day exposure data
Data_5day_old_2day_salinity_exposure <-
  Data_5day_old %>% 
  filter(`Larval salinity exposure duration` == "2 days")

Mean_bodydiameter_5day_old_2day_salinity_exposure <- Data_5day_old_2day_salinity_exposure %>% 
  group_by(`Larval salinity treatment`) %>% 
  summarize(
    mean_armdiam5 = mean(`Body diameter (arm tip-arm tip, um)`, na.rm = TRUE),
    std_error = sd(`Body diameter (arm tip-arm tip, um)`, na.rm = TRUE) / sqrt(n()),
    sample_size = n()
  )

print(Mean_bodydiameter_5day_old_2day_salinity_exposure)

#All 1 day old/4 day exposure data
Data_5day_old_4day_salinity_exposure <-
  Data_5day_old %>% 
  filter(`Larval salinity exposure duration` == "4 days")

mean_bodydiameter_5day_old_4day_salinity_exposure <- Data_5day_old_4day_salinity_exposure %>% 
  group_by(`Larval salinity treatment`) %>% 
  summarize(
    mean_armdiam6 = mean(`Body diameter (arm tip-arm tip, um)`, na.rm = TRUE),
    std_error = sd(`Body diameter (arm tip-arm tip, um)`, na.rm = TRUE) / sqrt(n()),
    sample_size = n()
  )

print(mean_bodydiameter_5day_old_4day_salinity_exposure)
```

#####Mean body diameter in each salinity, regardless of exposure time
```{r}
mean_bodydiameter_5day_old_sal <- Data_5day_old %>% 
  group_by(`Larval salinity treatment`) %>% 
  summarize(
    mean_armdiam8 = mean(`Body diameter (arm tip-arm tip, um)`, na.rm = TRUE),
    std_error = sd(`Body diameter (arm tip-arm tip, um)`, na.rm = TRUE) / sqrt(n()),
    sample_size = n()
  )

print(mean_bodydiameter_5day_old_sal)

#34-30
((669.4091-520.1389)/520.1389)*100

#34-22
((669.4091-420.7273)/420.7273)*100
```

#####Mean body diameter by exposure time, regardless of salinity
```{r}
#All 1 day old/2 day exposure data
Data_5day_old_2day_salinity_exposure <-
  Data_5day_old %>% 
  filter(`Larval salinity exposure duration` == "2 days")

mean_bodydiameter_5day_old_2day_salinity_exposure <- Data_5day_old_2day_salinity_exposure %>% 
  group_by(`Larval salinity exposure duration`) %>% 
  summarize(
    mean_armdiam1 = mean(`Body diameter (arm tip-arm tip, um)`, na.rm = T),
    std_error = sd(`Body diameter (arm tip-arm tip, um)`, na.rm = T) / sqrt(n()),
    sample_size = n()
  )

print(mean_bodydiameter_5day_old_2day_salinity_exposure)

#All 1 day old/4 day exposure data
Data_5day_old_4day_salinity_exposure <-
  Data_5day_old %>% 
  filter(`Larval salinity exposure duration` == "4 days")

mean_bodydiameter_5day_old_4day_salinity_exposure <- Data_5day_old_4day_salinity_exposure %>% 
  group_by(`Larval salinity exposure duration`) %>% 
  summarize(
    mean_armdiam1 = mean(`Body diameter (arm tip-arm tip, um)`, na.rm = TRUE),
    std_error = sd(`Body diameter (arm tip-arm tip, um)`, na.rm = TRUE) / sqrt(n()),
    sample_size = n()
  )

print(mean_bodydiameter_5day_old_4day_salinity_exposure)

#2d-4d
((582.0286-519.5577)/519.5577)*100
582.0286-519.5577
```

##Figure 2
```{r}
#Please note that some of the components of this figure as seen in the manuscript were added in photoshope software, after the base figures were generated in R (e.g. the figure legend)

# Calculate mean body diameter with standard error by salinity group
mean_bodydiameter_1d <- Data_1day_old %>%
  filter(!is.na(suppressWarnings(as.numeric(`Body diameter (arm tip-arm tip, um)`)))) %>%
  group_by(`Larval salinity treatment`, `Larval salinity exposure duration`) %>% 
  summarize(
    mean_bodydiam1 = mean(as.numeric(`Body diameter (arm tip-arm tip, um)`), na.rm = TRUE),
    std_error1 = sd(as.numeric(`Body diameter (arm tip-arm tip, um)`), na.rm = TRUE) / sqrt(n()),
    sample_size = n()
  )

# Print the calculated means and errors
print(mean_bodydiameter_1d)

# Create the plot
p1 <- ggplot(mean_bodydiameter_1d, aes(x = `Larval salinity treatment`, y = mean_bodydiam1, fill = `Larval salinity exposure duration`)) +
  geom_col(position = position_dodge(width = 0.8), width = 0.7) +
  geom_errorbar(aes(ymin = mean_bodydiam1 - std_error1, ymax = mean_bodydiam1 + std_error1),
                position = position_dodge(width = 0.8), width = 0.2) +
  scale_fill_manual(values = c("#d9d9d9", "#737373")) +
  ggtitle("", "1-day old juveniles") +
  labs(
    x = "Larval salinity treatment (â°)",
    y = "Mean diameter (mm, Â±se)",
    fill = "Salinity duration"
  ) +
  theme(
    plot.title = element_text(hjust = 0.0, vjust = 1, size = 12, family = "Times New Roman"),
    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.text = element_text(size = 12, family = "Times New Roman"),
    legend.title = element_text(size = 12, family = "Times New Roman"),
    legend.position = ""
  ) +
  scale_y_continuous(limits = c(0, 700), breaks = seq(0, 700, by = 100))

p1

# Calculate mean arm diameter with standard error by salinity group
mean_bodydiameter_5d <- Data_5day_old %>%
  filter(!is.na(suppressWarnings(as.numeric(`Body diameter (arm tip-arm tip, um)`)))) %>%  # Omit non-numeric cells
  group_by(`Larval salinity treatment`, `Larval salinity exposure duration`) %>% 
  summarize(
    mean_bodydiam2 = mean(as.numeric(`Body diameter (arm tip-arm tip, um)`), na.rm = TRUE),
    std_error1 = sd(as.numeric(`Body diameter (arm tip-arm tip, um)`), na.rm = TRUE) / sqrt(n()),
    sample_size = n()
  )

# Print the calculated means and errors
print(mean_bodydiameter_5d)

p2 <- ggplot(mean_bodydiameter_5d, aes(x = `Larval salinity treatment`, y = mean_bodydiam2, fill = `Larval salinity exposure duration`)) +
  geom_col(position = position_dodge(width = 0.8), width = 0.7) +
  geom_errorbar(aes(ymin = mean_bodydiam2 - std_error1, ymax = mean_bodydiam2 + std_error1),
                position = position_dodge(width = 0.8), width = 0.2) +
  scale_fill_manual(values = c("#d9d9d9", "#737373")) +
  ggtitle("", "5-day old juveniles") +
  labs(
    x = "Larval salinity treatment (â°)",
    y = "",
    fill = "Salinity duration"
  ) +
  theme(
    plot.title = element_text(hjust = 0.0, vjust = 1, size = 12, family = "Times New Roman"),
    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.text = element_text(size = 12, family = "Times New Roman"),
    legend.title = element_text(size = 12, family = "Times New Roman"),
    legend.position = ""
  ) +
  ylim(0, 700) +
  scale_y_continuous(breaks = seq(0, 700, by = 100))

p2

```

####Combined plots
```{r}
bodydiameter<-plot_grid(p1,p2, ncol = 2)
bodydiameter
```

#Influence of the larval salinity experience on juvenile morphology 
##Figure 3
###Load data
```{r}
library(dplyr)
Data_morphology <- read_excel("/Users/mcle5454/Desktop/For submission/Brach project 2/Post revision/Teratogenic effects of larval low salinity experience in development of the juvenile body in Acanthaster sp..xlsx", 
                           sheet = "Proportion normal", 
                           range = "A2:F41", 
                           col_names = c("Juvenile age",	"Larval salinity treatment",	"Larval salinity exposure duration",	"Normal category",	"prop value",	"Value (%)"))

# Transform relevant columns to factors
Data_morphology$`Larval salinity treatment` <- as.factor(Data_morphology$`Larval salinity treatment`)
Data_morphology$`Larval salinity exposure duration` <- as.factor(Data_morphology$`Larval salinity exposure duration`)
Data_morphology$`Normal category` <- as.factor(Data_morphology$`Normal category`)
```

###Graph
```{r}
#Stacked column graph (normal morphology, abnormal morphology and mortality)
##Please note that some of the components of this figure as seen in the manuscript were added in photoshope software, after the base figures were generated in R (e.g. the figure legend)

# Create the ggplot with cumulative SE
greyscale_colors <- c("#f0f0f0", "#d9d9d9", "#bdbdbd", "#969696", "#737373", "#525252", "#252525")

#2 d

Data_morphology_1day_old_2day_salinity_exposure<-
  Data_morphology %>% 
  filter(`Larval salinity exposure duration` =="2d", `Juvenile age` == "1 day old")

Data_morphology_1day_old_2day_salinity_exposure$`Normal category`<-factor(Data_morphology_1day_old_2day_salinity_exposure$`Normal category`,
                                                levels = c("Abnormal","Normal"))

Data_morphology_1day_old_2day_salinity_exposure$`Value (%)` <- as.numeric(Data_morphology_1day_old_2day_salinity_exposure$`Value (%)`)


p2.1 <- ggplot(Data_morphology_1day_old_2day_salinity_exposure, aes(x = `Larval salinity treatment`, y = `Value (%)`, fill = `Normal category`)) +
  geom_bar(stat = 'identity', position = "stack") +  # Stacking bars
  ggtitle("2 day larval salinity exposure","1 day old juveniles") +
  labs(x = "", y = "Proportion (%)", fill = NULL) +
  theme(
    plot.title = element_text(hjust = 0.0, vjust = 1, 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.position = ""
  ) +
  scale_fill_manual(values = c("#bdbdbd", "#737373")) +
  ylim(0, 100)  # Set y-axis limits

p2.1

#prop normal vs the original # of juveniles
Data_morphology_5day_old_2day_salinity_exposure<-
  Data_morphology %>% 
  filter(`Larval salinity exposure duration`=="2d", `Juvenile age` == "5 day old")

Data_morphology_5day_old_2day_salinity_exposure$`Normal category`<-factor(Data_morphology_5day_old_2day_salinity_exposure$`Normal category`,
                                                levels = c("Mortality","Abnormal","Normal"))

Data_morphology_5day_old_2day_salinity_exposure$`Value (%)` <- as.numeric(Data_morphology_5day_old_2day_salinity_exposure$`Value (%)`)


p2.2 <- ggplot(Data_morphology_5day_old_2day_salinity_exposure, aes(x = `Larval salinity treatment`, y = `Value (%)`, fill = `Normal category`)) +
  geom_bar(stat = 'identity', position = "stack") +  # Stacking bars
  ggtitle("","5 day old juveniles") +
  labs(x = "", y = "", fill = NULL) +
  theme(
    plot.title = element_text(hjust = 0.0, vjust = 1, 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.position = ""
  ) +
  scale_fill_manual(values = c("#f0f0f0", "#bdbdbd", "#737373")) +
  ylim(0, 100)  # Set y-axis limits

p2.2

#prop normal vs the remaining # of juveniles
Data_morphology_1day_old_4day_salinity_exposure<-
  Data_morphology %>% 
  filter(`Larval salinity exposure duration`=="4d", `Juvenile age` == "1 day old")

Data_morphology_1day_old_4day_salinity_exposure$`Normal category`<-factor(Data_morphology_1day_old_4day_salinity_exposure$`Normal category`,
                                                levels = c("Abnormal","Normal"))

Data_morphology_1day_old_4day_salinity_exposure$`Value (%)` <- as.numeric(Data_morphology_1day_old_4day_salinity_exposure$`Value (%)`)

p2.3 <- ggplot(Data_morphology_1day_old_4day_salinity_exposure, aes(x = `Larval salinity treatment`, y = `Value (%)`, fill = `Normal category`)) +
  geom_bar(stat = 'identity', position = "stack") +  # Stacking bars
  ggtitle("4 day larval salinity exposure","") +
  labs(x = "Larval salinity treatment (â°)", y = "Proportion (%)", fill = NULL) +
  theme(
    plot.title = element_text(hjust = 0.0, vjust = 1, 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.position = ""
  ) +
  scale_fill_manual(values = c("#bdbdbd", "#737373")) +
  ylim(0, 100)  # Set y-axis limits

p2.3

#prop normal vs the remaining # of juveniles
Data_morphology_5day_old_4day_salinity_exposure<-
  Data_morphology %>% 
  filter(`Larval salinity exposure duration`=="4d", `Juvenile age` == "5 day old")

Data_morphology_5day_old_4day_salinity_exposure$`Normal category`<-factor(Data_morphology_5day_old_4day_salinity_exposure$`Normal category`,
                                                levels = c("Mortality","Abnormal","Normal"))

Data_morphology_5day_old_4day_salinity_exposure$`Value (%)` <- as.numeric(Data_morphology_5day_old_4day_salinity_exposure$`Value (%)`)


p2.4 <- ggplot(Data_morphology_5day_old_4day_salinity_exposure, aes(x = `Larval salinity treatment`, y = `Value (%)`, fill = `Normal category`)) +
  geom_bar(stat = 'identity', position = "stack") +  # Stacking bars
  ggtitle("","") +
  labs(x = "Larval salinity treatment (â°)", y = "", fill = NULL) +
  theme(
    plot.title = element_text(hjust = 0.0, vjust = 1, 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.position = "",
  ) +
  scale_fill_manual(values = c("#f0f0f0", "#bdbdbd", "#737373")) +
  ylim(0, 100)  # Set y-axis limits

p2.4

```

###Combined plots
```{r}
stackedbps_norm_mort<-plot_grid(p2.1,p2.2,p2.3,p2.4, ncol = 2, nrow = 2)
stackedbps_norm_mort
```

##Logistic regressions and logistic regression plots 
###1 day old juveniles
###Load data 
```{r}
Data_morph_1dayold <- read_excel("/Users/mcle5454/Desktop/For submission/Brach project 2/Post revision/Teratogenic effects of larval low salinity experience in development of the juvenile body in Acanthaster sp..xlsx", 
                           sheet = "Juvenile morphology scores", 
                           range = "A2:M104", 
                           col_names = c("Ref #",	"Juvie Id",	"Larval salinity treatment",	"Larval salinity exposure duration",	"Replicate#",	"Juvenile#",	"Juvenile age (days)",	"Arm #",	"Normal (pentaradial)",	"Abnormal",	"4 arm",	"Fleshy/poor skeletal development",	"Swollen"))

# Transform relevant columns to factors
Data_morph_1dayold$sal <- as.factor(Data_morph_1dayold$`Larval salinity treatment`)
Data_morph_1dayold$Time <- as.factor(Data_morph_1dayold$`Larval salinity exposure duration`)
Data_morph_1dayold$vial <- as.factor(Data_morph_1dayold$`Replicate#`)#the larval vial replicate 
```

####GLMM_Morphology analyses
```{r}
glmer_control <- glmerControl(optimizer = "bobyqa")#(gold standard for use with repeated measures or convergence issues)
glmer_control <- glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 100000))

y_norm1<-Data_morph_1dayold$`Normal (pentaradial)`

Mod1<-glmer(y_norm1 ~ sal*Time+(1|vial), family = binomial(link = "logit"), data = Data_morph_1dayold, nAGQ = 25, control = glmer_control)
Aov1<-Anova(Mod1)
Aov1#not the best model fit

Mod2<-glmer(y_norm1 ~ sal+Time+(1|vial), family = binomial(link = "logit"), data = Data_morph1d, nAGQ = 25, control = glmer_control)
Aov2<-Anova(Mod2)
Aov2

emmip(Mod2, ~sal, CI=TRUE, type = "response") 

emmeans_Mod2<-emmeans(Mod2, ~sal, type = "response")#used for the probability curve

pairs <- pairs(emmeans_Mod2)
pairs

emmeans_Mod2<-emmeans(Mod2, ~sal, type = "response")#for curve
noperfect_separation_juv_2<-as.data.frame(emmeans_Mod2)

Probs_Juv_with_normalmorphology_1dayold <- ggplot(noperfect_separation_juv_2, aes(x = sal, y = prob, group = 1)) +  # Use linetype for differentiation
  geom_line(position = position_dodge(width = 0.25), size = 1, colour = "black") +  # Lines in black
  geom_point(position = position_dodge(width = 0.25), size = 3, colour = "black") +  # Points in black
  geom_errorbar(aes(ymin = asymp.LCL, ymax = asymp.UCL), 
                width = 0,        # Removes caps
                size = 0.6,      # Error bars slightly thicker
                alpha = 0.7,     # Keeps error bars 30% transparent
                position = position_dodge(width = 0.25), 
                colour = "black") +  # Error bars in black
  ylab("Prob(Juvenile with normal morphology)") +
  xlab("Larval salinity treatment (â°)") +
  ggtitle("1-day old juveniles") +
  theme(plot.title = element_text(hjust = 0.5, size = 12, family = "Times New Roman"),
        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"),
        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 = "",
        legend.text = element_text(size = 12, family = "Times New Roman"),
        legend.title = element_text(size = 12, family = "Times New Roman")) +
  ylim(0, 1)

Probs_Juv_with_normalmorphology_1dayold
```
#####Assumptions & data reliability checks
```{r}
#The following code tests the assumptions of the above analysis titled "LMM_1day_old (body diameter)":
##Residual diagnostics (to check for model misspecification (e.g., unaccounted nonlinearity or heteroscedasticity))
##Overdispersion (Tests whether residual variance exceeds that expected under the model)
##Zero-inflation (Checks whether there are more zeros than expected under the model distribution)
##Uniformity (KolmogorovâSmirnov test for global model fit & tests if simulated residuals are uniformly distributed)
##Random effect visualization (Visual check for normality and variation of random effects)
##Model discrimination and calibration (Checks calibration (agreement between predicted and observed responses)
##Model comparison (Likelihood ratio test to evaluate whether fixed effects improve model fit over random effects alone)

# Simulate residuals
simulation_output <- simulateResiduals(fittedModel = Mod1, n = 1000)

# Plot residuals
plot(simulation_output)  # Provides diagnostics for residuals (normality, dispersion, etc.)

# Test for overdispersion
testDispersion(simulation_output)

# Test for zero-inflation (if relevant)
testZeroInflation(simulation_output)

# Test for uniformity (global goodness-of-fit test)
testUniformity(simulation_output)

# Plot random effects
dotplot(ranef(Mod1, condVar = TRUE))  # Random effects with confidence intervals

# Create binned residual plots for each predictor
binnedplot(fitted(Mod1), residuals(Mod1, type = "response"), xlab = "Fitted Values", main = "Binned Residual Plot")

# Fit a null model without fixed effects
null_model <- glmer(y_norm1 ~ (1 | vial), family = binomial(link = "logit"), data = Data_morph_1dayold)

# Likelihood Ratio Test
anova(null_model, Mod1, test = "Chisq")
```
###5 day old juveniles
###Load data 
```{r}
Data_morph_5dayold <- read_excel("/Users/mcle5454/Desktop/For submission/Brach project 2/Post revision/Teratogenic effects of larval low salinity experience in development of the juvenile body in Acanthaster sp..xlsx", 
                           sheet = "Juvenile morphology scores", 
                           range = "A105:M207", 
                           col_names = c("Ref #",	"Juvie Id",	"Larval salinity treatment",	"Larval salinity exposure duration",	"Replicate#",	"Juvenile#",	"Juvenile age (days)",	"Arm #",	"Normal (pentaradial)",	"Abnormal",	"4 arm",	"Fleshy/poor skeletal development",	"Swollen"))

Data_morph_5dayold <- Data_morph_5dayold %>%
  filter(!is.na(suppressWarnings(as.numeric(`Normal (pentaradial)`)))) %>%
  mutate(`Normal (pentaradial)` = as.numeric(`Normal (pentaradial)`))

# Transform relevant columns to factors
Data_morph_5dayold$sal <- as.factor(Data_morph_5dayold$`Larval salinity treatment`)
Data_morph_5dayold$Time <- as.factor(Data_morph_5dayold$`Larval salinity exposure duration`)
Data_morph_5dayold$vial <- as.factor(Data_morph_5dayold$`Replicate#`)#larval vial replicate
```

####GLMM_Morphology analyses
```{r}
glmer_control <- glmerControl(optimizer = "bobyqa")#(gold standard for use with repeated measures or convergence issues)
glmer_control <- glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 100000))

y_norm2<-Data_morph_5dayold$`Normal (pentaradial)`

Mod3<-glmer(y_norm2 ~ sal*Time+(1|vial), family = binomial(link = "logit"), data = Data_morph_5dayold, nAGQ = 25, control = glmer_control)#not the best model fit
Aov3<-Anova(Mod3)
Aov3

Mod4<-glmer(y_norm2 ~ sal+Time+(1|vial), family = binomial(link = "logit"), data = Data_morph_5dayold, nAGQ = 25, control = glmer_control)
Aov4<-Anova(Mod4)
Aov4

emmip(Mod4, ~sal, CI=TRUE, type = "response") 
emmeans_Mod4<-emmeans(Mod4, ~sal, type = "response")

pairs <- pairs(emmeans_Mod4)
pairs

emmeans_Mod4<-emmeans(Mod4, ~sal, type = "response")#for curve
noperfect_separation_juv_4<-as.data.frame(emmeans_Mod4)

Probs_Juv_with_normalmorphology_5dayold <- ggplot(noperfect_separation_juv_4, aes(x = sal, y = prob, group = 1)) +  # Use linetype for differentiation
  geom_line(position = position_dodge(width = 0.25), size = 1, colour = "black") +  # Lines in black
  geom_point(position = position_dodge(width = 0.25), size = 3, colour = "black") +  # Points in black
  geom_errorbar(aes(ymin = asymp.LCL, ymax = asymp.UCL), 
                width = 0,        # Removes caps
                size = 0.6,      # Error bars slightly thicker
                alpha = 0.7,     # Keeps error bars 30% transparent
                position = position_dodge(width = 0.25), 
                colour = "black") +  # Error bars in black
  ylab("") +
  xlab("Larval salinity treatment (â°)") +
  ggtitle("5-day old juveniles") +
  theme(plot.title = element_text(hjust = 0.5, size = 12, family = "Times New Roman"),
        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"),
        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 = "",
        legend.text = element_text(size = 12, family = "Times New Roman"),
        legend.title = element_text(size = 12, family = "Times New Roman")) +
  ylim(0, 1)

Probs_Juv_with_normalmorphology_5dayold
```
#####Assumptions & data reliability checks
```{r}
#The following code tests the assumptions of the above analysis titled "LMM_1day_old (body diameter)":
##Residual diagnostics (to check for model misspecification (e.g., unaccounted nonlinearity or heteroscedasticity))
##Overdispersion (Tests whether residual variance exceeds that expected under the model)
##Zero-inflation (Checks whether there are more zeros than expected under the model distribution)
##Uniformity (KolmogorovâSmirnov test for global model fit & tests if simulated residuals are uniformly distributed)
##Random effect visualization (Visual check for normality and variation of random effects)
##Model discrimination and calibration (Checks calibration (agreement between predicted and observed responses)
##Model comparison (Likelihood ratio test to evaluate whether fixed effects improve model fit over random effects alone)

# Simulate residuals
simulation_output <- simulateResiduals(fittedModel = Mod2, n = 1000)

# Plot residuals
plot(simulation_output)  # Provides diagnostics for residuals (normality, dispersion, etc.)

# Test for overdispersion
testDispersion(simulation_output)

# Test for zero-inflation (if relevant)
testZeroInflation(simulation_output)

# Test for uniformity (global goodness-of-fit test)
testUniformity(simulation_output)

# Plot random effects
dotplot(ranef(Mod2, condVar = TRUE))  # Random effects with confidence intervals

# Create binned residual plots for each predictor
binnedplot(fitted(Mod2), residuals(Mod2, type = "response"), xlab = "Fitted Values", main = "Binned Residual Plot")

# Fit a null model without fixed effects
null_model <- glmer(y_norm2 ~ (1 | vial), family = binomial(link = "logit"), data = Data_morph_5dayold)

# Likelihood Ratio Test
anova(null_model, Mod2, test = "Chisq")
```


###(Figure 4)
#####Combined plots of GLMM logistic regressions above
```{r}
##Binomial logistic regression curves for the probability of a juvenile having a normal morphology with experience to 2 or 4 days of exposure to salinity levels: 34 to 22 ppt. Please note that some of the components of this figure as seen in the manuscript were added in photoshope software, after the base figures were generated in R (e.g. the figure legend)

Juv_morphology<-plot_grid(Probs_Juv_with_normalmorphology_1dayold, Probs_Juv_with_normalmorphology_5dayold, ncol = 2)
Juv_morphology
```

###Arm#
####Load data
```{r}
library(dplyr)
ArmData <- read_excel("/Users/mcle5454/Desktop/For submission/Brach project 2/Post revision/Teratogenic effects of larval low salinity experience in development of the juvenile body in Acanthaster sp..xlsx", 
                           sheet = "Proportional arm count", 
                           range = "A2:F49", 
                           col_names = c("Juvenile age",	"Larval salinity treatment",	"Larval salinity exposure duration",	"Arm number",	"Proportion",	"Percentage (%)"))

# Transform relevant columns to factors
ArmData$`Larval sal` <- as.factor(ArmData$`Larval salinity treatment`)
ArmData$`Days exposed` <- as.factor(ArmData$`Larval salinity exposure duration`)
ArmData$`Arm number` <- as.factor(ArmData$`Arm number`)
```

####Figure 5
```{r}
# Create the ggplot with cumulative SE
greyscale_colors <- c("#f0f0f0", "#d9d9d9", "#bdbdbd", "#969696", "#737373", "#525252", "#252525")

#2 d

ArmData_1day_old_2day_salinity_exposure<-
  ArmData %>% 
  filter(`Days exposed`=="2d", `Juvenile age` == "1 day old")

ArmData_1day_old_2day_salinity_exposure$`Arm number`<-as.factor(ArmData_1day_old_2day_salinity_exposure$`Arm number`)

df3 <- ArmData_1day_old_2day_salinity_exposure %>%
  group_by(`Larval salinity treatment`, `Arm number`) %>%
  summarise(mean3 = mean(`Percentage (%)`),
            .groups = 'drop') %>%
  mutate(`Larval salinity treatment` = factor(`Larval salinity treatment`, levels = c(22, 25, 30, 34)),
        `Arm number` = factor(`Arm number`))

p1 <- ggplot(df3, aes(x = `Larval salinity treatment`, y = mean3, fill = `Arm number`)) +
  geom_bar(stat = 'identity') + 
  ggtitle("2 day larval salinity exposure","1 day old juveniles") +
  labs(x = "", y = "Proportion (%)", fill = NULL) +
  theme(plot.title = element_text(hjust = 0.0, vjust = 1, 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.position = "") +  # Adjust legend position to the right
  scale_fill_manual(values = c("#f0f0f0", "#bdbdbd", "#525252")) +
  ylim(0, 100) +
  scale_y_continuous(limits = c(0, 100),breaks = c(0, 20, 40, 60, 80, 100))

p1



ArmData_5day_old_2day_salinity_exposure<-
  ArmData %>% 
  filter(`Days exposed`=="2d", `Juvenile age` == "5 day old")

ArmData_1day_old_2day_salinity_exposure$`Arm number`<-as.factor(ArmData_1day_old_2day_salinity_exposure$`Arm number`)

df5 <- ArmData_1day_old_2day_salinity_exposure %>%
  group_by(`Larval salinity treatment`, `Arm number`) %>%
  summarise(mean5 = mean(`Percentage (%)`),
            .groups = 'drop') %>%
  mutate(`Larval salinity treatment` = factor(`Larval salinity treatment`, levels = c(22, 25, 30, 34)),
        `Arm number` = factor(`Arm number`))

p2 <- ggplot(df5, aes(x = `Larval salinity treatment`, y = mean5, fill = `Arm number`)) +
  geom_bar(stat = 'identity') + 
  ggtitle("","5 day old juveniles") +
  labs(x = "", y = "", fill = NULL) +
  theme(plot.title = element_text(hjust = 0.0, vjust = 1, size = 12, family = "Times New Roman"),
        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.position = "") +  # Adjust legend position to the right
  scale_fill_manual(values = c("#f0f0f0", "#bdbdbd", "#525252")) +
  ylim(0, 100) +
  scale_y_continuous(limits = c(0, 100),breaks = c(0, 20, 40, 60, 80, 100))

p2



ArmData_1day_old_4day_salinity_exposure<-
  ArmData %>% 
  filter(`Days exposed`=="4d", `Juvenile age` == "1 day old")

ArmData_1day_old_4day_salinity_exposure$`Arm number`<-as.factor(ArmData_1day_old_4day_salinity_exposure$`Arm number`)

df6 <- ArmData_1day_old_4day_salinity_exposure %>%
  group_by(`Larval salinity treatment`, `Arm number`) %>%
  summarise(mean6 = mean(`Percentage (%)`),
            .groups = 'drop') %>%
  mutate(`Larval salinity treatment` = factor(`Larval salinity treatment`, levels = c(22, 25, 30, 34)),
        `Arm number` = factor(`Arm number`))

p3 <- ggplot(df6, aes(x = `Larval salinity treatment`, y = mean6, fill = `Arm number`)) +
  geom_bar(stat = 'identity') + 
  ggtitle("4 day larval salinity exposure","") +
  labs(x = "Larval salinity treatment (â°)", y = "Proportion (%)", fill = NULL) +
  theme(plot.title = element_text(hjust = 0.0, vjust = 1, 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.position = "") +  # Adjust legend position to the right
  scale_fill_manual(values = c("#f0f0f0", "#bdbdbd", "#525252")) +
  ylim(0, 100) +
  scale_y_continuous(limits = c(0, 100),breaks = c(0, 20, 40, 60, 80, 100))

p3



ArmData_5day_old_4day_salinity_exposure<-
  ArmData %>% 
  filter(`Days exposed`=="4d", `Juvenile age` == "5 day old")

df7 <- ArmData_5day_old_4day_salinity_exposure %>%
  group_by(`Larval salinity treatment`, `Arm number`) %>%
  summarise(mean7 = mean(`Percentage (%)`),
            .groups = 'drop') %>%
  mutate(`Larval salinity treatment` = factor(`Larval salinity treatment`, levels = c(22, 25, 30, 34)),
        `Arm number` = factor(`Arm number`))

p4 <- ggplot(df7, aes(x = `Larval salinity treatment`, y = mean7, fill = `Arm number`)) +
  geom_bar(stat = 'identity') + 
  ggtitle("","") +
  labs(x = "Larval salinity treatment (â°)", y = "", fill = NULL) +
  theme(plot.title = element_text(hjust = 0.0, vjust = 1, size = 12, family = "Times New Roman"),
        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.position = "") +  # Adjust legend position to the right
  scale_fill_manual(values = c("#f0f0f0", "#bdbdbd", "#525252")) +
  ylim(0, 100) +
  scale_y_continuous(limits = c(0, 100),breaks = c(0, 20, 40, 60, 80, 100))

p4
```

#####Combined plots
```{r}
##Stacked column graph (no arms, 4 arms, and 5 arms). Please note that some of the components of this figure as seen in the manuscript were added in photoshope software, after the base figures were generated in R (e.g. the figure legend)

Arm_no<-plot_grid(p1,p2,p3,p4, ncol = 2, nrow = 2)
Arm_no
```

##Analyses
###Load Data (1 day old) 
```{r}
Data_arm_1dayold <- read_excel("/Users/mcle5454/Desktop/For submission/Brach project 2/Post revision/Teratogenic effects of larval low salinity experience in development of the juvenile body in Acanthaster sp..xlsx", 
                           sheet = "Arm#", 
                           range = "A2:K104", 
                           col_names = c("Ref #",	"Juvie Id",	"Larval salinity treatment",	"Larval salinity exposure duration",	"Replicate#",	"Juvenile#",	"Juvenile age (days)",	"Body diameter (arm tip-arm tip, um)",	"5 arms",	"4 arms",	"0 arms"))

# Transform relevant columns to factors
Data_arm_1dayold$sal <- as.factor(Data_arm_1dayold$`Larval salinity treatment`)
Data_arm_1dayold$Time <- as.factor(Data_arm_1dayold$`Larval salinity exposure duration`)
Data_arm_1dayold$vial <- as.factor(Data_arm_1dayold$`Replicate#`)
```

######5Arm analysis_day 1
```{r}
glmer_control <- glmerControl(optimizer = "bobyqa")#(gold standard for use with repeated measures or convergence issues)
glmer_control <- glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 100000))

arm5<-Data_arm_1dayold$`5 arms`

#run the rptd measures GLMM
Mod7<-glmer(arm5 ~ sal*Time+(1|vial), family = binomial(link = "logit"), data = Data_arm_1dayold, nAGQ = 25, control = glmer_control)
Aov7<-Anova(Mod7)
Aov7

Mod8<-glmer(arm5 ~ sal+Time+(1|vial), family = binomial(link = "logit"), data = Data_arm_1dayold, nAGQ = 25, control = glmer_control)
Aov8<-Anova(Mod8)
Aov8


emmeans_Mod8<-emmeans(Mod8, ~sal, type = "response")

pairs <- pairs(emmeans_Mod8)
pairs
```

#####Assumptions & data reliability checks
```{r}
#The following code tests the assumptions of the above analysis titled "LMM_1day_old (body diameter)":
##Residual diagnostics (to check for model misspecification (e.g., unaccounted nonlinearity or heteroscedasticity))
##Overdispersion (Tests whether residual variance exceeds that expected under the model)
##Zero-inflation (Checks whether there are more zeros than expected under the model distribution)
##Uniformity (KolmogorovâSmirnov test for global model fit & tests if simulated residuals are uniformly distributed)
##Random effect visualization (Visual check for normality and variation of random effects)
##Model discrimination and calibration (Checks calibration (agreement between predicted and observed responses)
##Model comparison (Likelihood ratio test to evaluate whether fixed effects improve model fit over random effects alone)

# Simulate residuals
simulation_output <- simulateResiduals(fittedModel = Mod8, n = 1000)

# Plot residuals
plot(simulation_output)  # Provides diagnostics for residuals (normality, dispersion, etc.)

# Test for overdispersion
testDispersion(simulation_output)

# Test for zero-inflation (if relevant)
testZeroInflation(simulation_output)

# Test for uniformity (global goodness-of-fit test)
testUniformity(simulation_output)

# Fit a GLM (without random effects) to approximate VIF
fixed_effects_model <- glm(arm5 ~ sal + Time, family = binomial(link = "logit"), data = Data_arm_1dayold)
vif(fixed_effects_model)

# Plot random effects
dotplot(ranef(Mod8, condVar = TRUE))  # Random effects with confidence intervals

# Create binned residual plots for each predictor
binnedplot(fitted(Mod8), residuals(Mod8, type = "response"), xlab = "Fitted Values", main = "Binned Residual Plot")

# Fit a null model without fixed effects
null_model <- glmer(arm5 ~ (1 | vial), family = binomial(link = "logit"), data = Data_arm_1dayold)

# Likelihood Ratio Test
anova(null_model, Mod8, test = "Chisq")

# Predicted probabilities
Data_arm_1dayold$predicted <- predict(Mod8, type = "response")

# Create a plot
ggplot(Data_arm_1dayold, aes(x = predicted, y = arm5)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "loess", se = FALSE) +
  labs(x = "Predicted Probability", y = "Observed Response")

```
### Load Data (5 day old)
```{r}
Data_arm_5dayold <- read_excel("/Users/mcle5454/Desktop/For submission/Brach project 2/Post revision/Teratogenic effects of larval low salinity experience in development of the juvenile body in Acanthaster sp..xlsx", 
                           sheet = "Arm#", 
                           range = "A105:K191", 
                           col_names = c("Ref #",	"Juvie Id",	"Larval salinity treatment",	"Larval salinity exposure duration",	"Replicate#",	"Juvenile#",	"Juvenile age (days)",	"Body diameter (arm tip-arm tip, um)",	"5 arms",	"4 arms",	"0 arms"))

# Transform relevant columns to factors
Data_arm_5dayold$sal <- as.factor(Data_arm_5dayold$`Larval salinity treatment`)
Data_arm_5dayold$Time <- as.factor(Data_arm_5dayold$`Larval salinity exposure duration`)
Data_arm_5dayold$vial <- as.factor(Data_arm_5dayold$`Replicate#`)
```

######5Arm analysis_day 5
```{r}
glmer_control <- glmerControl(optimizer = "bobyqa")#(gold standard for use with repeated measures or convergence issues)
glmer_control <- glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 100000))

arm5<-Data_arm_5dayold$`5 arms`

#run the rptd measures GLMM
Mod9<-glmer(arm5 ~ sal*Time+(1|vial), family = binomial(link = "logit"), data = Data_arm_5dayold, nAGQ = 25, control = glmer_control)
Aov9<-Anova(Mod9)
Aov9

Mod10<-glmer(arm5 ~ sal+Time+(1|vial), family = binomial(link = "logit"), data = Data_arm_5dayold, nAGQ = 25, control = glmer_control)
Aov10<-Anova(Mod10)
Aov10

emmeans_Mod10<-emmeans(Mod10, ~sal, type = "response")#for curve
pairs <- pairs(emmeans_Mod10)
pairs
```

#####Assumptions & data reliability checks
```{r}
#The following code tests the assumptions of the above analysis titled "LMM_1day_old (body diameter)":
##Residual diagnostics (to check for model misspecification (e.g., unaccounted nonlinearity or heteroscedasticity))
##Overdispersion (Tests whether residual variance exceeds that expected under the model)
##Zero-inflation (Checks whether there are more zeros than expected under the model distribution)
##Uniformity (KolmogorovâSmirnov test for global model fit & tests if simulated residuals are uniformly distributed)
##Random effect visualization (Visual check for normality and variation of random effects)
##Model discrimination and calibration (Checks calibration (agreement between predicted and observed responses)
##Model comparison (Likelihood ratio test to evaluate whether fixed effects improve model fit over random effects alone)

# Simulate residuals
simulation_output <- simulateResiduals(fittedModel = Mod10, n = 1000)

# Plot residuals
plot(simulation_output)  # Provides diagnostics for residuals (normality, dispersion, etc.)

# Test for overdispersion
testDispersion(simulation_output)

# Test for zero-inflation (if relevant)
testZeroInflation(simulation_output)

# Test for uniformity (global goodness-of-fit test)
testUniformity(simulation_output)

# Fit a GLM (without random effects) to approximate VIF
fixed_effects_model <- glm(arm5 ~ sal + Time, family = binomial(link = "logit"), data = Data_arm_5dayold)
vif(fixed_effects_model)

# Plot random effects
dotplot(ranef(Mod10, condVar = TRUE))  # Random effects with confidence intervals

# Create binned residual plots for each predictor
binnedplot(fitted(Mod10), residuals(Mod8, type = "response"), xlab = "Fitted Values", main = "Binned Residual Plot")

# Fit a null model without fixed effects
null_model <- glmer(arm5 ~ (1 | vial), family = binomial(link = "logit"), data = Data_arm_5dayold)

# Likelihood Ratio Test
anova(null_model, Mod10, test = "Chisq")

# Predicted probabilities
Data_arm_5dayold$predicted <- predict(Mod10, type = "response")

# Create a plot
ggplot(Data_arm_5dayold, aes(x = predicted, y = arm5)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "loess", se = FALSE) +
  labs(x = "Predicted Probability", y = "Observed Response")

```

