---
title: "Chapter 1 _ MHW impacts on juveniles"
author: "J Donelson and Y Yasutake"
date: "Created Jan 2024, Updated Jan 2026"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)

knitr::opts_knit$set(root.dir = "/Users/jc492396/Library/CloudStorage/OneDrive-JamesCookUniversity/JCU PhD/Chapters/Chapter 1_ Juv growth/Stats/FINAL/Data")
```

## R Markdown
```{r Package loading, message=FALSE}
library(ggplot2)
library(tidyverse)
library(lme4)
library(lmerTest)
library(emmeans)
library(car)
library(performance)
library(readxl)
```

```{r , include=FALSE}
getwd()
```

# Data loading
```{r}
Egg <- read_excel(
  path = "Embryonic exposure to a marine heatwave restricts juvenile growth in a coral reef damselfish.xlsx", 
                     sheet = "Egg size", 
                     trim_ws = TRUE)

Egg$Treatment = factor(Egg$Treatment)
Egg$Clutch = factor(Egg$Clutch)
Egg$Parental.tank = factor(Egg$Parental.tank)
Egg$Male = factor(Egg$Male)
Egg$Female = factor(Egg$Female)
Egg$Sample = factor(Egg$Sample)
Egg$Egg.size = as.numeric(Egg$Egg.size)

str(Egg)
```

```{r}
MHW <- read_excel(
  path = "Embryonic exposure to a marine heatwave restricts juvenile growth in a coral reef damselfish.xlsx", 
                     sheet = "Day 0-60", 
                     trim_ws = TRUE)

MHW$Clutch = factor(MHW$Clutch)
MHW$Treatment = factor(MHW$Treatment)
MHW$Day = factor(MHW$Day)
MHW$Male = factor(MHW$Male)
MHW$Female = factor(MHW$Female)
MHW$Fish.rep = factor(MHW$Fish.rep)
MHW$Parental.tank = factor(MHW$Parental.tank)
MHW$Tank = factor(MHW$Tank)

MHW$SL = as.numeric(MHW$SL)
MHW$W  = as.numeric(MHW$W)
MHW$YA  = as.numeric(MHW$YA)
MHW$FemaleSL  = as.numeric(MHW$FemaleSL)
MHW$FemaleW  = as.numeric(MHW$FemaleW)
MHW$FemaleFK  = as.numeric(MHW$FemaleFK)
MHW$MaleSL  = as.numeric(MHW$MaleSL)
MHW$MaleW  = as.numeric(MHW$MaleW)
MHW$MaleFK  = as.numeric(MHW$MaleFK)

MHW$LogFemaleSL <- log(MHW$FemaleSL)
MHW$LogMaleSL <- log(MHW$MaleSL)
MHW$LogFemaleFK <- log(MHW$FemaleFK)
MHW$LogMaleFK <- log(MHW$MaleFK)
MHW$LogSL <- log(MHW$SL)
MHW$LogW <- log(MHW$W)

str(MHW)
```


```{r}
MHWembryo <- MHW %>%
  filter(Treatment != "Early exposure",
         Treatment != "Late exposure",
         Day != 0)

str(MHWembryo)
```
```{r}
MHWday15_30 <- MHW %>%
  filter(Day %in% c("15", "30"),
         Treatment != "Embryo",
         !Parental.tank %in% c("31", "144", "216")) # removing pairs that are not present in all treatments

str(MHWday15_30)
```

```{r}
MHWday45_60 <- MHW %>%
  filter(Day %in% c("45", "60"),
         Treatment != "Embryo",
         !Parental.tank %in% c("31", "144", "216")) # removing pairs that are not present in all treatments

str(MHWday45_60)
```

```{r}
Hatching = MHW %>%
  filter(Day %in% c("0"))

str(Hatching)
```

```{r}
EarlyExposure <- MHW %>%
  filter(Day %in% c("15", "30", "45", "60"),
         Treatment %in% c("Control", "Early exposure"),
         !Parental.tank %in% c("31", "144", "216"))

str(EarlyExposure)
```

# Overview
For each trait, information are given as follows:  
1. Random effect structures  
2. Covariate exploration  
3. Final model  
4. Model performance  
5. Model outputs (summary statistics, anova, emmeans, and pairwise comparisons where appropriate)  
6. Figure  

# Effects of MHWs during embryonic development  
# Egg size (Control vs Embryo)  
## 1. Random effect structure  
```{r}
lmer.egg = lmer(Egg.size ~ Treatment + (1|Parental.tank/Clutch), data = Egg)
```
Clutch nested within Parental tank as 2 clutches from Pair 32 and 66 are used (1 each in Control and Embryo).

## 2. Covariate exploration
```{r}
lm = lm(Egg.size ~ FemaleSL, data=Egg)
summary(lm)

lm = lm(Egg.size ~ FemaleFK, data=Egg)
summary(lm)
```

Stronger correlation with maternal FK. Testing model fit.  

```{r}
lmer.egg = lmer(Egg.size ~ Treatment + (1|Parental.tank/Clutch), data = Egg)
lmer.egga = lmer(Egg.size ~ Treatment + FemaleFK + (1|Parental.tank/Clutch), data = Egg)

AIC(lmer.egg, lmer.egga)
```

Choosing a (with FemaleFK)

## 3. Final model 
```{r}
lmer.egg = lmer(Egg.size ~ Treatment + FemaleFK + (1|Parental.tank/Clutch), data = Egg)
```

## 4. Performance check
```{r}
performance::check_model(lmer.egg, check="homogeneity") 
performance::check_model(lmer.egg, check="outliers") 
performance::check_model(lmer.egg, check="qq", detrend = FALSE) 
performance::check_model(lmer.egg, check="normality") 
performance::check_model(lmer.egg, check="linearity") 
performance::check_model(lmer.egg, check="pp_check")
hist(residuals(lmer.egg), col="darkgray")
shapiro.test(residuals(lmer.egg))
qqline(resid(lmer.egg))
outlierTest(lmer.egg)
```

\#92 identified as outlier but the egg is normal (e.g., not squished). Kept. 

## 5. Outputs
### Summary statistics
```{r}
summary(lmer.egg)
```

### anova 
```{r}
Anova(lmer.egg, type = 3)
```
### Emmeans
```{r}
EmmEgg = (emmeans(lmer.egg, ~ Treatment) %>% as.data.frame) 
EmmEgg
```

## 6. Figure
```{r}
fig.Egg <- ggplot(EmmEgg, 
                  aes(x = Treatment, y = emmean, fill = Treatment)) + 
  geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL),
                position = position_dodge(width = 0.6),
                width = 0.1, colour = "black", linewidth = 0.3) + 
  geom_point(shape = 22, size = 8, colour = "black", stroke = 0.3,
             position = position_dodge(width = 0.6)) +
  theme_classic() +
  labs(y = "Egg size (mm²)") + 
  scale_fill_manual(values = c("Control" = "#0072B2", 
                               "Embryo" = "#B2182B")) +
  scale_y_continuous(
  limits = c(4, 5.5),
  breaks = c(4, 4.5, 5, 5.5),
  labels = c(4, 4.5, 5, 5.5)) +
  theme(
  text = element_text(family = "Helvetica"),
  axis.title.x = element_blank(),
  axis.text.x  = element_blank(),
  axis.title.y = element_text(size = 14, face = "bold", colour = "black"),
  axis.text.y  = element_text(size = 12, colour = "black"),
  legend.position = "none"
)

fig.Egg
```

# Juvenile morphometrics at hatching
## Standard length at hatching
## 1. Random effect structures
```{r}
lmer.HatchSL = lmer(SL ~ Treatment + (1|Clutch), data = Hatching)
lmer.HatchSLa = lmer(SL ~ Treatment  + (1|Parental.tank/Clutch), data = Hatching)

AIC(lmer.HatchSL, lmer.HatchSLa)
```

Choosing "a" to account for cases where clutches from the same pair is used.  

## 2. Covariate exploration
```{r}
lm = lm(SL ~ FemaleSL, data=Hatching)
summary(lm)

lm = lm(SL ~ FemaleFK, data=Hatching)
summary(lm)

lm = lm(SL ~ MaleSL, data=Hatching)
summary(lm)

lm = lm(SL ~ MaleFK, data=Hatching)
summary(lm)
```
Strongest correlation with Female SL. Testing model fit.  

```{r}
lmer.HatchSLa = lmer(SL ~ Treatment  + (1|Parental.tank/Clutch), data = Hatching)
lmer.HatchSLb = lmer(SL ~ Treatment  + FemaleSL + (1|Parental.tank/Clutch), data = Hatching)

AIC(lmer.HatchSLa, lmer.HatchSLb)
```
Choosing a (without covariate).

## 3. Final model 
```{r}
lmer.HatchSL = lmer(SL ~ Treatment  + (1|Parental.tank/Clutch), data = Hatching)
```

## 4. Model performance
```{r}
performance::check_model(lmer.HatchSL, check="homogeneity") 
performance::check_model(lmer.HatchSL, check="outliers") 
performance::check_model(lmer.HatchSL, check="qq", detrend = FALSE) 
performance::check_model(lmer.HatchSL, check="normality") 
performance::check_model(lmer.HatchSL, check="linearity") 
performance::check_model(lmer.HatchSL, check="pp_check", re_formula=NA)
hist(residuals(lmer.HatchSL), col="darkgray")
shapiro.test(residuals(lmer.HatchSL))
qqline(resid(lmer.HatchSL))
outlierTest(lmer.HatchSL)
```

\#320 identified as an outlier but retained as it appears within the expected range for that clutch.   

## 5. Output
### Summary statistics
```{r}
summary(lmer.HatchSL)
```

### anova
```{r}
Anova(lmer.HatchSL, type = 3)
```

### Emmean
```{r}
EmmHatchSL = (emmeans(lmer.HatchSL, ~ Treatment) %>% as.data.frame) 

EmmHatchSL
```
## 6. Figure
```{r}
fig.HatchSL <- ggplot(EmmHatchSL, 
                     aes(x = Treatment, y = emmean, fill = Treatment)) + 
  geom_errorbar(
    aes(ymin = lower.CL, ymax = upper.CL),
    position = position_dodge(width = 0.6),
    width = 0.1, colour = "black", linewidth = 0.3
  ) + 
  geom_point(
    shape = 22, size = 8, colour = "black", stroke = 0.3,
    position = position_dodge(width = 0.6)
  ) +
  theme_classic() +
  labs(
    x = NULL,
    y = "Standard legnth (cm)"
    ) +
  scale_fill_manual(values = c("Control" = "#0072B2", 
                               "Embryo" = "#B2182B")) +    
  scale_y_continuous(
  limits = c(0.5, 0.6),
  breaks = c(0.5, 0.55, 0.60),
  labels = c("0.50", "0.550", "0.60")) +
  theme(
    text = element_text(family = "Helvetica"),
    axis.title = element_text(size = 14, face = "bold", colour = "black"),
    axis.text = element_text(size = 12, colour = "black"),
    axis.text.x  = element_blank(),
    legend.position = "none"
  )

fig.HatchSL
```

## Juvenile weight at hatching
## 1. Random effect structure
```{r}
lmer.HatchW = lmer(W ~ Treatment + (1|Clutch), data = Hatching)
lmer.HatchWa = lmer(W ~ Treatment  + (1|Parental.tank/Clutch), data = Hatching)

AIC(lmer.HatchW, lmer.HatchWa)
```
As per above, choosing "a" to account for same clutches from the same pair.  

## 2. Covariates
```{r}
lm = lm(W ~ FemaleSL, data=Hatching)
summary(lm)

lm = lm(W ~ FemaleFK, data=Hatching)
summary(lm)

lm = lm(W ~ MaleSL, data=Hatching)
summary(lm)

lm = lm(W ~ MaleFK, data=Hatching)
summary(lm)
```

Strongest correlation with Female SL. Testing model fit. 

```{r}
lmer.HatchWa = lmer(W ~ Treatment  + (1|Parental.tank/Clutch), data = Hatching)
lmer.HatchWb = lmer(W ~ Treatment  + FemaleSL + (1|Parental.tank/Clutch), data = Hatching) # fails to converge
```

b fails to converge. Choosing a (without covariate).

## 3. Final model 
```{r}
lmer.HatchW = lmer(W ~ Treatment  + (1|Parental.tank/Clutch), data = Hatching)
```

## 4. Model performance
```{r}
performance::check_model(lmer.HatchW, check="homogeneity") 
performance::check_model(lmer.HatchW, check="outliers") 
performance::check_model(lmer.HatchW, check="qq", detrend = FALSE) 
performance::check_model(lmer.HatchW, check="normality") 
performance::check_model(lmer.HatchW, check="linearity") 
performance::check_model(lmer.HatchW, check="pp_check", re_formula=NA)
hist(residuals(lmer.HatchW), col="darkgray")
shapiro.test(residuals(lmer.HatchW))
qqline(resid(lmer.HatchW))
outlierTest(lmer.HatchW)
```

\#89 identified as an outlier but retained as it appears within expected range of that clutch.  

## 5. Output
### Summary statistics
```{r}
summary(lmer.HatchW)
```

### anova
```{r}
Anova(lmer.HatchW, type = 3)
```

### Emmean
```{r}
EmmHatchW = (emmeans(lmer.HatchW, ~ Treatment) %>% as.data.frame) 

EmmHatchW
```
## 6. Figure
```{r}
fig.HatchW <- ggplot(EmmHatchW, aes(x = Treatment, y = emmean, fill = Treatment)) + 
  geom_errorbar(
    aes(ymin = lower.CL, ymax = upper.CL),
    position = position_dodge(width = 0.6),
    width = 0.1, colour = "black", linewidth = 0.3
  ) + 
  geom_point(
    shape = 22, size = 8, colour = "black", stroke = 0.3,
    position = position_dodge(width = 0.6)
  ) +
  theme_classic() +
  labs(
    x = NULL, 
    y = "Weight (g)"
  ) +
  scale_fill_manual(values = c("Control" = "#0072B2", "Embryo" = "#B2182B")) +
  scale_y_continuous(
  limits = c(0.002, 0.004),
  breaks = c(0.002, 0.003, 0.004),
  labels = c(0.002, 0.003, 0.004)) +
  theme(
    text = element_text(family = "Helvetica"),
    axis.title.x = element_text(size = 14, face = "bold", colour = "black", margin = margin(t = 8)),
    axis.title.y = element_text(size = 14, face = "bold", colour = "black", margin = margin(r = 8)),
    axis.text = element_text(size = 12, colour = "black"),
    axis.text.x  = element_blank(),
    legend.position = "none"
  )

fig.HatchW
```

# Condition (W for given SL) at hatching 
## 1. Random effect structure
```{r}
lmer.HatchCond = lmer(W ~ Treatment + SL + (1|Clutch), data = Hatching)
lmer.HatchConda = lmer(W ~ Treatment  + SL + (1|Parental.tank/Clutch), data = Hatching)

AIC(lmer.HatchCond, lmer.HatchConda)
```

Choosing "a" to account for cases where clutches from the same pair is used.  

## 2. Covariate exploration
Testing Female SL to follow hatching W.

```{r}
lmer.HatchConda = lmer(W ~ Treatment  + SL + (1|Parental.tank/Clutch), data = Hatching)
lmer.HatchCondb = lmer(W ~ Treatment  + SL + FemaleSL + (1|Parental.tank/Clutch), data = Hatching)

AIC(lmer.HatchConda, lmer.HatchCondb)
```

Choosing a (without covariate). 

## 3. Final model 
```{r}
lmer.HatchCond = lmer(W ~ Treatment  + SL + (1|Parental.tank/Clutch), data = Hatching)
```

## 4. Model performance
```{r}
performance::check_model(lmer.HatchCond, check="homogeneity") 
performance::check_model(lmer.HatchCond, check="outliers") 
performance::check_model(lmer.HatchCond, check="qq", detrend = FALSE) 
performance::check_model(lmer.HatchCond, check="normality") 
performance::check_model(lmer.HatchCond, check="linearity") 
performance::check_model(lmer.HatchCond, check="pp_check", re_formula=NA)
hist(residuals(lmer.HatchCond), col="darkgray")
shapiro.test(residuals(lmer.HatchCond))
qqline(resid(lmer.HatchCond))
outlierTest(lmer.HatchCond)
```

\#155 identified as an outlier but retained as it appears within the range for that clutch.  

## 5. Output
### Summary statistics 
```{r}
summary(lmer.HatchCond)
```

### anova
```{r}
Anova(lmer.HatchCond, type = 3)
```

### Emmean
```{r}
EmmHatchCond = (emmeans(lmer.HatchCond, ~ Treatment) %>% as.data.frame) 

EmmHatchCond
```
## 6. Figure
```{r}
fig.HatchCond <- ggplot(EmmHatchCond, 
                     aes(x = Treatment, y = emmean, fill = Treatment)) + 
  geom_errorbar(
    aes(ymin = lower.CL, ymax = upper.CL),
    position = position_dodge(width = 0.6),
    width = 0.1, colour = "black", linewidth = 0.3
  ) + 
  geom_point(
    shape = 22, size = 8, colour = "black", stroke = 0.3,
    position = position_dodge(width = 0.6)
  ) +
  theme_classic() +
  labs(
    y = "Condition (g)"
  ) +
  scale_y_continuous(
  limits = c(0.002, 0.004),
  breaks = c(0.002, 0.003, 0.004),
  labels = c(0.002, 0.003, 0.004)) +
  scale_fill_manual(values = c("Control" = "#0072B2", 
                               "Embryo" = "#B2182B")) + 
  theme(
    text = element_text(family = "Helvetica"),
    axis.title.x = element_blank(),
    axis.title = element_text(size = 14, face = "bold", colour = "black"),
    axis.text = element_text(size = 12, colour = "black"),
    legend.position = "none"
  )

fig.HatchCond
```

## Yolk area at hatching
## 1.  Random effect structure
```{r}
lmer.HatchYA = lmer(YA ~ Treatment + (1|Clutch), data = Hatching)
lmer.HatchYAa = lmer(YA ~ Treatment  + (1|Parental.tank/Clutch), data = Hatching)

AIC(lmer.HatchYA, lmer.HatchYAa)
```

Choosing "a" to account for cases where clutches from the same pair is used.  
 

## 2. Covariates
```{r}
lm = lm(YA ~ FemaleSL, data=Hatching)
summary(lm)

lm = lm(YA ~ FemaleFK, data=Hatching)
summary(lm)

lm = lm(YA ~ MaleSL, data=Hatching)
summary(lm)

lm = lm(YA ~ MaleFK, data=Hatching)
summary(lm)
```

Strongest correlation with Female SL. Testing model fit. 

```{r}
lmer.HatchYAa = lmer(YA ~ Treatment  + (1|Parental.tank/Clutch), data = Hatching)
lmer.HatchYAb = lmer(YA ~ Treatment  + FemaleSL + (1|Parental.tank/Clutch), data = Hatching)

AIC(lmer.HatchYAa, lmer.HatchYAb)
```

Choosing a (without covariate).  

## 3. Final model 
```{r}
lmer.HatchYA = lmer(YA ~ Treatment  + (1|Parental.tank/Clutch), data = Hatching)
```

## 4. Model performance
```{r}
performance::check_model(lmer.HatchYA, check="homogeneity") 
performance::check_model(lmer.HatchYA, check="outliers") 
performance::check_model(lmer.HatchYA, check="qq", detrend = FALSE) 
performance::check_model(lmer.HatchYA, check="normality") 
performance::check_model(lmer.HatchYA, check="linearity") 
performance::check_model(lmer.HatchYA, check="pp_check", re_formula=NA)
hist(residuals(lmer.HatchYA), col="darkgray")
shapiro.test(residuals(lmer.HatchYA))
qqline(resid(lmer.HatchYA))
outlierTest(lmer.HatchYA)
```

\#345 identified as an outlier but retained as it appears within the expected range of the clutch.  

## 5. Model outputs
### Summary statistics
```{r}
summary(lmer.HatchYA)
```

### anova
```{r}
Anova(lmer.HatchYA, type = 3)
```

### Emmean
```{r}
EmmHatchYA = (emmeans(lmer.HatchYA, ~ Treatment) %>% as.data.frame) 

EmmHatchYA
```
## 6. Figure 
```{r}
fig.HatchYA <- ggplot(EmmHatchYA, 
                     aes(x = Treatment, y = emmean, fill = Treatment)) + 
  geom_errorbar(
    aes(ymin = lower.CL, ymax = upper.CL),
    position = position_dodge(width = 0.6),
    width = 0.1, colour = "black", linewidth = 0.3
  ) + 
  geom_point(
    shape = 22, size = 8, colour = "black", stroke = 0.3,
    position = position_dodge(width = 0.6)
  ) +
  theme_classic() +
  labs(
    x = NULL,
    y = "Yolk area (cm²)"
    ) +
  scale_fill_manual(values = c("Control" = "#0072B2", 
                               "Embryo" = "#B2182B")) +    
  scale_y_continuous(
  limits = c(0.01, 0.02),
  breaks = c(0.01, 0.015,0.02),
  labels = c(0.01, 0.015,0.02)) +
  theme(
    text = element_text(family = "Helvetica"),
    axis.title = element_text(size = 14, face = "bold", colour = "black"),
    axis.text = element_text(size = 12, colour = "black"),
    legend.position = "none"
    )

fig.HatchYA
```

Below are effects of MHWs during embryonic development on post-hatching growth. 

# Growth in standard length
## 1. Random effect structures

```{r}
lmer.SLa = lmer(LogSL ~ Treatment * Day + (1|Tank) + (1|Parental.tank/Clutch), data = MHWembryo)
lmer.SLb = lmer(LogSL ~ Treatment * Day + (1|Parental.tank/Clutch/Tank), data = MHWembryo)

AIC(lmer.SLa, lmer.SLb)
```

lmer.SL and lmer.SLb essentially a same model. Opting for the nested model. 

## 2. Covariate exploration
```{r}
lm = lm(SL ~ FemaleSL, data=MHWembryo)
summary(lm)

lm = lm(SL ~ FemaleFK, data=MHWembryo)
summary(lm)

lm = lm(SL ~ MaleSL, data=MHWembryo)
summary(lm)

lm = lm(SL ~ MaleFK, data=MHWembryo)
summary(lm)
```

No parental trait with significant correlation. 

## 3. Final model 
```{r}
lmer.SL = lmer(LogSL ~ Treatment * Day + (1|Parental.tank/Clutch/Tank), data = MHWembryo)
```

## 4. Model performance
```{r}
performance::check_model(lmer.SL, check="homogeneity") 
performance::check_model(lmer.SL, check="outliers") 
performance::check_model(lmer.SL, check="qq", detrend = FALSE) 
performance::check_model(lmer.SL, check="normality") 
performance::check_model(lmer.SL, check="linearity") 
performance::check_model(lmer.SL, check="pp_check", re_formula=NA)
hist(residuals(lmer.SL), col="darkgray")
shapiro.test(residuals(lmer.SL))
qqline(resid(lmer.SL))
outlierTest(lmer.SL)
```

\#522 identified as outlier but retained as it appears within the expected range for that clutch.

## 5. Model outouts 
### Summary statistics
```{r}
summary(lmer.SL)
```

### anova
```{r}
Anova(lmer.SL, type = 3)
```
### Emmean
```{r}
EmmSL.Log = (emmeans(lmer.SL, ~ Treatment * Day) %>% as.data.frame) 

# Back-transforming
EmmSL <- EmmSL.Log %>%
  mutate(
    emmean = exp(emmean),
    lower.CL = exp(lower.CL),
    upper.CL = exp(upper.CL)
  )

EmmSL
```

### Pairwise comparison
#### Day
```{r}
emm_day <- emmeans(lmer.SL, ~ Day)

pairs_day <- pairs(emm_day, adjust = "tukey") |> summary(infer = TRUE)
pairs_day
```

## 2. Figure
```{r}
EmmSL$Day <- as.numeric(as.character(EmmSL$Day))

fig.Embryo.SL <- ggplot(EmmSL,
       aes(x = Day, y = emmean, colour = Treatment, fill = Treatment)) +
  geom_ribbon(aes(ymin = lower.CL, ymax = upper.CL),
              alpha = 0.15, colour = NA) +
  geom_line(linewidth = 1) +
  geom_point(size = 3, stroke = 0.5, colour = "black",
             aes(fill = Treatment), shape = 22) +

  scale_fill_manual(values = c("Control" = "#0072B2",
                               "Embryo"  = "#B2182B")) +
  scale_colour_manual(values = c("Control" = "#0072B2",
                                 "Embryo"  = "#B2182B")) +
  scale_x_continuous(breaks = c(15, 30, 45, 60),
                     labels = c(15, 30, 45, 60)) +
  scale_y_continuous(
  limits = c(0.5, 2.5),
  breaks = c(0.5, 1, 1.5, 2, 2.5),
  labels = c(0.5, 1, 1.5, 2, 2.5)) +
  theme_classic() +
  labs(
    y = "Standard length (cm)",
    x = NULL
  ) +
  theme(
    text = element_text(family = "Helvetica"),
    axis.title = element_text(size = 14, face = "bold", colour = "black"),
    axis.text = element_text(size = 12, colour = "black"),
    axis.text.x  = element_blank(),
    legend.position = "none"
    )

fig.Embryo.SL
```

# Growth in Weight
## 1. Random effect structure
```{r}
lmer.Wa = lmer(LogW ~ Treatment * Day + (1|Tank) + (1|Parental.tank/Clutch), data = MHWembryo)
lmer.Wb = lmer(LogW ~ Treatment * Day + (1|Parental.tank/Clutch/Tank), data = MHWembryo)

AIC(lmer.Wa, lmer.Wb)
```
Choosing b. 

## 2. Covariate exploration
```{r}
lm = lm(W ~ FemaleSL, data=MHWembryo)
summary(lm)

lm = lm(W ~ FemaleFK, data=MHWembryo)
summary(lm)

lm = lm(W ~ MaleSL, data=MHWembryo)
summary(lm)

lm = lm(W ~ MaleFK, data=MHWembryo)
summary(lm)
```

No parental trait with significant correlation. 

## 3. Final model
```{r}
lmer.W = lmer(LogW ~ Treatment * Day + (1|Parental.tank/Clutch/Tank), data = MHWembryo)
```

## 4. Model performance
```{r}
performance::check_model(lmer.W, check="homogeneity") 
performance::check_model(lmer.W, check="outliers") 
performance::check_model(lmer.W, check="qq", detrend = FALSE) 
performance::check_model(lmer.W, check="normality") 
performance::check_model(lmer.W, check="linearity") 
performance::check_model(lmer.W, check="pp_check")
hist(residuals(lmer.W), col="darkgray")
shapiro.test(residuals(lmer.W))
qqline(resid(lmer.W))
outlierTest(lmer.W)
```
\#552 identified as an outlier but retained as it appears within the range for that clutch.

## 5. Model output
### Summary statistics 
```{r}
summary(lmer.W)
```

### anova
```{r}
Anova(lmer.W, type = 3)
```
## Pairwise comparison
### Day
```{r}
emm_day <- emmeans(lmer.W, ~ Day)

pairs_day <- pairs(emm_day, adjust = "tukey") |> summary(infer = TRUE)
pairs_day
```
### Emmean
```{r}
EmmW.Log <- emmeans(lmer.W, ~ Treatment * Day) %>% 
  as.data.frame()

EmmW <- EmmW.Log %>%
  mutate(
    emmean   = exp(emmean),
    lower.CL = exp(lower.CL),
    upper.CL = exp(upper.CL)
  )

EmmW
```

## 6. Figure
```{r}
EmmW$Day <- as.numeric(as.character(EmmW$Day))

fig.Embryo.W <- ggplot(EmmW,
       aes(x = Day, y = emmean, colour = Treatment, fill = Treatment)) +
  geom_ribbon(aes(ymin = lower.CL, ymax = upper.CL),
              alpha = 0.15, colour = NA) +
  geom_line(linewidth = 1) +
  geom_point(size = 3, stroke = 0.5, colour = "black",
             aes(fill = Treatment), shape = 22) +

  scale_fill_manual(values = c("Control" = "#0072B2",
                               "Embryo"  = "#B2182B")) +
  scale_colour_manual(values = c("Control" = "#0072B2",
                                 "Embryo"  = "#B2182B")) +
  scale_x_continuous(breaks = c(15, 30, 45, 60),
                     labels = c(15, 30, 45, 60)) +
    scale_y_continuous(
  limits = c(0, 0.45),
  breaks = c(0, 0.1, 0.2, 0.3, 0.4),
  labels = c(0, 0.1, 0.2, 0.3, 0.4)) +

  theme_classic() +
  labs(
    y = "Weight (g)",
    x = NULL
  ) +
  theme(
    text = element_text(family = "Helvetica"),
    axis.title = element_text(size = 14, face = "bold", colour = "black"),
    axis.text = element_text(size = 12, colour = "black"),
    axis.text.x  = element_blank(),
    legend.position = "none"
    )

fig.Embryo.W
```



# Growth in Condition (W for given SL)
```{r}
MHWembryo <- MHWembryo %>%
  group_by(Day) %>%
  mutate(LogSL.centered = LogSL - mean(LogSL, na.rm = TRUE))
```

## 1. Random effect structure
```{r}
lmer.Conda = lmer(LogW ~ Treatment * Day + LogSL.centered + (1|Tank) + (1|Parental.tank/Clutch), data = MHWembryo) # fails to converge
lmer.Condb = lmer(LogW ~ Treatment * Day + LogSL.centered + (1|Parental.tank/Clutch/Tank), data = MHWembryo) # fails to converge
lmer.Condc = lmer(LogW ~ Treatment * Day + LogSL.centered + (1|Clutch/Tank), data = MHWembryo) 
```
a and b fails to converge. Choosing c. 

## 2. Covariate explotation
No parental trait with significant correlation for W. 

## 3. Final model 
```{r}
lmer.Cond = lmer(LogW ~ Treatment * Day + LogSL.centered + (1|Clutch/Tank), data = MHWembryo) 
```

## 4. Model performance
```{r}
performance::check_model(lmer.Cond, check="homogeneity") 
performance::check_model(lmer.Cond, check="outliers") 
performance::check_model(lmer.Cond, check="qq", detrend = FALSE) 
performance::check_model(lmer.Cond, check="normality") 
performance::check_model(lmer.Cond, check="linearity") 
performance::check_model(lmer.Cond, check="pp_check")
hist(residuals(lmer.Cond), col="darkgray")
shapiro.test(residuals(lmer.Cond))
qqline(resid(lmer.Cond))
outlierTest(lmer.Cond)
```

All outliers appeared to be within the expected range for a particular clutch. Kept. 

## 5. Model outputs
### Summary statistics 
```{r}
summary(lmer.Cond)
```

### anova
```{r}
Anova(lmer.Cond, type = 3)
```
### Emmean
```{r}
EmmCond.Log <- emmeans(lmer.Cond, ~ Treatment * Day) |> as.data.frame()

EmmCond <- EmmCond.Log |>
  mutate(
    emmean   = exp(emmean),
    lower.CL = exp(lower.CL),
    upper.CL = exp(upper.CL)
  )

EmmCond
```

## Pariwise
### Day 
```{r}
emm_Cond <- emmeans(lmer.Cond, ~ Day)

pairs_Cond <- pairs(emm_Cond, adjust = "tukey") |> summary(infer = TRUE)
pairs_Cond
```

### Treatment x Day
```{r}
emm.Cond <- emmeans(lmer.Cond, ~ Treatment | Day)

pairs.Cond <- pairs(emm.Cond, adjust = "tukey") %>%
  summary(infer = TRUE)

pairs.Cond
```

## 6. Figure
```{r}
EmmCond$Day <- as.numeric(as.character(EmmCond$Day))

fig.Embryo.Cond <- ggplot(EmmCond,
       aes(x = Day, y = emmean, colour = Treatment, fill = Treatment)) +
  geom_ribbon(aes(ymin = lower.CL, ymax = upper.CL),
              alpha = 0.15, colour = NA) +
  geom_line(linewidth = 1) +
  geom_point(size = 3, stroke = 0.5, colour = "black",
             aes(fill = Treatment), shape = 22) +

  scale_fill_manual(values = c("Control" = "#0072B2",
                               "Embryo"  = "#B2182B")) +
  scale_colour_manual(values = c("Control" = "#0072B2",
                                 "Embryo"  = "#B2182B")) +
  scale_x_continuous(breaks = c(15, 30, 45, 60),
                     labels = c(15, 30, 45, 60)) +
  scale_y_continuous(
  limits = c(0, 0.4),
  breaks = c(0, 0.1, 0.2, 0.3, 0.4),
  labels = c(0, 0.1, 0.2, 0.3, 0.4)) +
  theme_classic() +
  labs(
    y = "Condition (g)",
    x = NULL
  ) +
  theme(
    text = element_text(family = "Helvetica"),
    axis.title = element_text(size = 14, face = "bold", colour = "black"),
    axis.text = element_text(size = 12, colour = "black"),
    legend.position = "none"
    )

fig.Embryo.Cond
```

# Effects of post-hatching MHWs on juvenile growth
# Grwoth in SL (Day 15 and 30)

## 1. Random effect structures
Each clutch came from unique parents unlike impacts of embryonic MHWs above. Parental.tank are not included in the models below for this reason. 

```{r}
lmer.SL15a = lmer(LogSL ~ Treatment * Day + (1|Tank) + (1|Clutch), data = MHWday15_30)
lmer.SL15b = lmer(LogSL ~ Treatment * Day + (1|Clutch/Tank), data = MHWday15_30)

AIC(lmer.SL15a, lmer.SL15b)
```

Choosing b. 

## 2. Covariate exploration
```{r}
lm = lm(SL ~ FemaleSL, data=MHWday15_30)
summary(lm)

lm = lm(SL ~ FemaleFK, data=MHWday15_30)
summary(lm)

lm = lm(SL ~ MaleSL, data=MHWday15_30)
summary(lm)

lm = lm(SL ~ MaleFK, data=MHWday15_30)
summary(lm)
```

Strongest correlation with Female SL. Testing model fit. 

```{r}
lmer.SL15b = lmer(LogSL ~ Treatment * Day + (1|Clutch/Tank), data = MHWday15_30)
lmer.SL15c = lmer(LogSL ~ Treatment * Day + LogFemaleSL + (1|Clutch/Tank), data = MHWday15_30)

AIC(lmer.SL15b, lmer.SL15c)
```
Choosing c (with LogFemaleSL). 

## 3. Final model
```{r}
lmer.SL15 = lmer(LogSL ~ Treatment * Day + LogFemaleSL + (1|Clutch/Tank), data = MHWday15_30)
```

## 4. Model performance
```{r}
performance::check_model(lmer.SL15, check="homogeneity") 
performance::check_model(lmer.SL15, check="outliers") 
performance::check_model(lmer.SL15, check="qq", detrend = FALSE) 
performance::check_model(lmer.SL15, check="normality") 
performance::check_model(lmer.SL15, check="linearity") 
performance::check_model(lmer.SL15, check="pp_check")
hist(residuals(lmer.SL15), col="darkgray")
shapiro.test(residuals(lmer.SL15))
qqline(resid(lmer.SL15))
outlierTest(lmer.SL15)
```

All outliers checked and removed if they were unusual for a given clutch.

## 5. Model outputs
### Summary statistics
```{r}
summary(lmer.SL15)
```

### anova
```{r}
Anova(lmer.SL15, type = 3)
```

### Emmean
```{r}
EmmSL.15_Log = (emmeans(lmer.SL15, ~ Treatment * Day) %>% as.data.frame) 

# Back-transforming
EmmSL.15_30 <- EmmSL.15_Log %>%
  mutate(
    emmean = exp(emmean),
    lower.CL = exp(lower.CL),
    upper.CL = exp(upper.CL)
  )

EmmSL.15_30
```

### Pairwise
```{r}
emm_day <- emmeans(lmer.SL15, ~ Day)

pairs_day <- pairs(emm_day, adjust = "tukey") |> summary(infer = TRUE)
pairs_day
```

## 6. Figure 
```{r}
fig.SL.15_30 <- ggplot(EmmSL.15_30, 
                        aes(x = Day, y = emmean, fill = Treatment)) + 
  geom_errorbar(
    aes(ymin = lower.CL, ymax = upper.CL),
    position = position_dodge(width = 0.6),     
    width = 0.1, colour = "black", linewidth = 0.3
  ) + 
  geom_point(
  aes(size = interaction(Treatment, Day)),  
  shape = 22, colour = "black", stroke = 0.3,
  position = position_dodge(width = 0.6)
) +
scale_size_manual(
  values = c(
    "Control.15" = 6,
    "Early exposure.15"  = 6,
    "Control.30" = 6,
    "Early exposure.30"  = 6
  ),
  guide = "none"  
) +
  theme_classic() +
  labs(
    y = "Standard length (cm)"
  ) +
  scale_fill_manual(values = c("Control" = "#0072B2", 
                               "Early exposure" = "#F2C200"
                               )) +    
  scale_y_continuous(
  limits = c(0.5, 1.5),
  breaks = c(0.5, 1, 1.5),
  labels = c(0.5, 1, 1.5)) +
  theme(
    text = element_text(family = "Helvetica"),
    axis.title.x = element_blank(),
    axis.title = element_text(size = 14, face = "bold", colour = "black"),
    axis.text = element_text(size = 12, colour = "black"),
    axis.text.x = element_blank(),
    legend.position = "none"
  )

fig.SL.15_30
```

# Growth in weight (Day 15 and 30)
## 1. Random effect structures
```{r}
lmer.W15a = lmer(LogW ~ Treatment * Day + (1|Tank) + (1|Clutch), data = MHWday15_30)
lmer.W15b = lmer(LogW ~ Treatment * Day + (1|Clutch/Tank) , data = MHWday15_30)

AIC(lmer.W15a, lmer.W15b)
```

Choosing b. 

## 2. Covariate exploration
```{r}
lm = lm(W ~ FemaleSL, data=MHWday15_30)
summary(lm)

lm = lm(W ~ FemaleFK, data=MHWday15_30)
summary(lm)

lm = lm(W ~ MaleSL, data=MHWday15_30)
summary(lm)

lm = lm(W ~ MaleFK, data=MHWday15_30)
summary(lm)
```

Strongest correlation with Female SL. Testing model fit.

```{r}
lmer.W15b = lmer(LogW ~ Treatment * Day + (1|Clutch/Tank) , data = MHWday15_30)
lmer.W15c = lmer(LogW ~ Treatment * Day + LogFemaleSL + (1|Clutch/Tank) , data = MHWday15_30)

AIC(lmer.W15b, lmer.W15c)
```

Choosing c (with LogFemaleSL). 

## 3. Final model
```{r}
lmer.W15 = lmer(LogW ~ Treatment * Day + LogFemaleSL + (1|Clutch/Tank) , data = MHWday15_30)
```

## 4. Model performance
```{r}
performance::check_model(lmer.W15, check="homogeneity") 
performance::check_model(lmer.W15, check="outliers") 
performance::check_model(lmer.W15, check="qq", detrend = FALSE) 
performance::check_model(lmer.W15, check="normality") 
performance::check_model(lmer.W15, check="linearity") 
performance::check_model(lmer.W15, check="pp_check")
hist(residuals(lmer.W15), col="darkgray")
shapiro.test(residuals(lmer.W15))
qqline(resid(lmer.W15))
outlierTest(lmer.W15)
```
\#143 identified as an outlier but retained as it appears within the expected range for that clutch. 

## 5. Model outputs
### Summary statistics
```{r}
summary(lmer.W15)
```

### anova
```{r}
Anova(lmer.W15, type = 3)
```
### Emmean
```{r}
EmmW.15_Log = emmeans(lmer.W15, ~ Treatment * Day) %>% as.data.frame() 

# Back-transforming
EmmW.15_30 <- EmmW.15_Log %>%
  mutate(
    emmean = exp(emmean),
    lower.CL = exp(lower.CL),
    upper.CL = exp(upper.CL)
  )

EmmW.15_30
```

### Pairwise
```{r}
emm_day <- emmeans(lmer.W15, ~ Day)

pairs_day <- pairs(emm_day, adjust = "tukey") |> summary(infer = TRUE)
pairs_day
```

## 6. Figure
```{r}
fig.W.15_30 <- ggplot(EmmW.15_30, 
                        aes(x = Day, y = emmean, fill = Treatment)) + 
  geom_errorbar(
    aes(ymin = lower.CL, ymax = upper.CL),
    position = position_dodge(width = 0.6),     
    width = 0.1, colour = "black", linewidth = 0.3
  ) + 
  geom_point(
  aes(size = interaction(Treatment, Day)),  
  shape = 22, colour = "black", stroke = 0.3,
  position = position_dodge(width = 0.6)
) +
scale_size_manual(
  values = c(
    "Control.15" = 6,
    "Early exposure.15"  = 6,
    "Control.30" = 6,
    "Early exposure.30"  = 6
  ),
  guide = "none"  
) +
  theme_classic() +
  labs(
    y = "Weight (g)" 
  ) +
  scale_fill_manual(values = c("Control" = "#0072B2", 
                               "Early exposure" = "#F2C200"
                               )) +    
  scale_y_continuous(
  limits = c(0, 0.1),
  breaks = c(0, 0.05, 0.1),
  labels = c(0, 0.05, 0.1)) +
  theme(
    text = element_text(family = "Helvetica"),
    axis.title.x = element_blank(),
    axis.title = element_text(size = 14, face = "bold", colour = "black"),
    axis.text = element_text(size = 12, colour = "black"),
    axis.text.x = element_blank(),
    legend.position = "none"
    )

fig.W.15_30
```

# Growth in Condition (Day 15 and 30)
## 1. Random effect structure
```{r}
MHWday15_30 <- MHWday15_30 %>%
  group_by(Day) %>%
  mutate(LogSL.centered = LogSL - mean(LogSL, na.rm = TRUE))
```

```{r}
lmer.Cond15a = lmer(LogW ~ Treatment * Day + LogSL.centered  + (1|Tank) + (1|Clutch), data = MHWday15_30)
lmer.Cond15b = lmer(LogW ~ Treatment * Day + LogSL.centered  + (1|Clutch/Tank), data = MHWday15_30)

AIC(lmer.Cond15a, lmer.Cond15b)
```

Choosing b. 

## 2. Covariate exploration
Testing model fit with Female SL to follow the model for W above. 

```{r}
lmer.Cond15b = lmer(LogW ~ Treatment * Day + LogSL.centered  + (1|Clutch/Tank), data = MHWday15_30)
lmer.Cond15c = lmer(LogW ~ Treatment * Day + LogSL.centered + LogFemaleSL + (1|Clutch/Tank), data = MHWday15_30)

AIC(lmer.Cond15b, lmer.Cond15c)
```

Choosing b (without covariate). 

## 3. Final model 
```{r}
lmer.Cond15 = lmer(LogW ~ Treatment * Day + LogSL.centered  + (1|Clutch/Tank), data = MHWday15_30)
```

## 4. Model performance
```{r}
performance::check_model(lmer.Cond15, check="homogeneity") 
performance::check_model(lmer.Cond15, check="outliers") 
performance::check_model(lmer.Cond15, check="qq", detrend = FALSE) 
performance::check_model(lmer.Cond15, check="normality") 
performance::check_model(lmer.Cond15, check="linearity") 
performance::check_model(lmer.Cond15, check="pp_check")
hist(residuals(lmer.Cond15), col="darkgray")
shapiro.test(residuals(lmer.Cond15))
qqline(resid(lmer.Cond15))
outlierTest(lmer.Cond15)
```

\#339 identified as an outlier but retained as it appears within the expected range of the clutch. 

## 5. Model output
### Summary statistics
```{r}
summary(lmer.Cond15)
```

### anova
```{r}
Anova(lmer.Cond15, type = 3)
```

### Emmean
```{r}
EmmCond.15_Log = (emmeans(lmer.Cond15, ~ Treatment * Day) %>% as.data.frame) 

# Back-transforming
EmmCond.15_30 <- EmmCond.15_Log %>%
  mutate(
    emmean = exp(emmean),
    lower.CL = exp(lower.CL),
    upper.CL = exp(upper.CL)
  )

EmmCond.15_30
```

### Pairwise
#### Day
```{r}
emm_day <- emmeans(lmer.Cond15, ~ Day)

pairs_day <- pairs(emm_day, adjust = "tukey") |> summary(infer = TRUE)
pairs_day
```

## 6. Figure 
```{r}
fig.Cond.15_30 <- ggplot(EmmCond.15_30, 
                        aes(x = Day, y = emmean, fill = Treatment)) + 
  geom_errorbar(
    aes(ymin = lower.CL, ymax = upper.CL),
    position = position_dodge(width = 0.6),     
    width = 0.1, colour = "black", linewidth = 0.3
  ) + 
  geom_point(
  aes(size = interaction(Treatment, Day)),  
  shape = 22, colour = "black", stroke = 0.3,
  position = position_dodge(width = 0.6)
) +
scale_size_manual(
  values = c(
    "Control.15" = 4.5,
    "Early exposure.15"  = 4.5,
    "Control.30" = 6,
    "Early exposure.30"  = 6
  ),
  guide = "none"  
) +
  theme_classic() +
  labs(
    y = "Condition (g)"
  ) +
  scale_fill_manual(values = c("Control" = "#0072B2", 
                               "Early exposure" = "#F2C200"
                               )) +   
  scale_y_continuous(
  limits = c(0.02, 0.08),
  breaks = c(0.02, 0.04, 0.06, 0.08),
  labels = c(0.02, 0.04, 0.06, 0.08)) +
  theme(
    text = element_text(family = "Helvetica"),
    axis.title.x = element_blank(),
    axis.title = element_text(size = 14, face = "bold", colour = "black"),
    axis.text = element_text(size = 12, colour = "black"),
    legend.position = "none"
    )

fig.Cond.15_30
```

# Growth in SL (Day 45 and 60)
## 1. Random effect structures
```{r}
lmer.SL45a = lmer(SL ~ Treatment * Day + (1|Tank) + (1|Clutch), data = MHWday45_60)
lmer.SL45b = lmer(SL ~ Treatment * Day  + (1|Clutch/Tank), data = MHWday45_60)

AIC(lmer.SL45a, lmer.SL45b)
```

Choosing b (without covariate). 

## 2. Covariate exploration
```{r}
lm = lm(SL ~ FemaleSL, data=MHWday45_60)
summary(lm)

lm = lm(SL ~ FemaleFK, data=MHWday45_60)
summary(lm)

lm = lm(SL ~ MaleSL, data=MHWday45_60)
summary(lm)

lm = lm(SL ~ MaleFK, data=MHWday45_60)
summary(lm)
```

Correlation strongest with Male SL Testing model fit.

```{r}
lmer.SL45b = lmer(SL ~ Treatment * Day  + (1|Clutch/Tank), data = MHWday45_60)
lmer.SL45c = lmer(SL ~ Treatment * Day  + MaleSL + (1|Clutch/Tank), data = MHWday45_60)

AIC(lmer.SL45b, lmer.SL45c)
```
Choosing b (without covariate). 

## 3. Final model
```{r}
lmer.SL45 = lmer(SL ~ Treatment * Day + (1|Clutch/Tank), data = MHWday45_60)
```

## 4. Performance check
```{r}
performance::check_model(lmer.SL45, check="homogeneity") 
performance::check_model(lmer.SL45, check="outliers") 
performance::check_model(lmer.SL45, check="qq", detrend = FALSE) 
performance::check_model(lmer.SL45, check="normality") 
performance::check_model(lmer.SL45, check="linearity") 
performance::check_model(lmer.SL45, check="pp_check")
hist(residuals(lmer.SL45), col="darkgray")
qqline(resid(lmer.SL45))
outlierTest(lmer.SL45)
```
\#80 identified as an outlier but retained as it appears within expected range for that clutch. 

## 5. Outputs
### Summary
```{r}
summary(lmer.SL45)
```

### anova
```{r}
Anova(lmer.SL45, type = 3)
```
### Emmeans
```{r}
EmmSL.45_60 = (emmeans(lmer.SL45, ~ Treatment * Day) %>% as.data.frame) 
EmmSL.45_60
```
### Pairwise
#### Day 
```{r}
emm_day <- emmeans(lmer.SL45, ~ Day)

pairs_day <- pairs(emm_day, adjust = "tukey") |> summary(infer = TRUE)
pairs_day
```

## 6. Figure
```{r}
fig.SL.45_60 <- ggplot(EmmSL.45_60, 
                        aes(x = Day, y = emmean, fill = Treatment)) + 
  geom_errorbar(
    aes(ymin = lower.CL, ymax = upper.CL),
    position = position_dodge(width = 0.6),     
    width = 0.1, colour = "black", linewidth = 0.3
  ) + 
  geom_point(
  aes(size = interaction(Treatment, Day)),  
  shape = 22, colour = "black", stroke = 0.3,
  position = position_dodge(width = 0.6)
) +
scale_size_manual(
  values = c(
    "Control.45" = 6,
    "Early exposure.45"  = 6,
    "Late exposure.45"  = 6,
    "Control.60" = 6,
    "Early exposure.60"  = 6,
    "Late exposure.60"  = 6
  ),
  guide = "none"  
) +
  theme_classic() +
  labs(
    y = "Standard length (cm)" 
  ) +
  scale_fill_manual(values = c("Control" = "#0072B2", 
                               "Early exposure" = "#F2C200",
                               "Late exposure" = "#D55E00")) +   
  scale_y_continuous(
  limits = c(1.5, 2.5),
  breaks = c(1.5, 2, 2.5),
  labels = c(1.5, 2, 2.5)) +
  theme(
    text = element_text(family = "Helvetica"),
    axis.title.x = element_blank(),
    axis.title = element_text(size = 14, face = "bold", colour = "black"),
    axis.text = element_text(size = 12, colour = "black"),
    axis.text.x  = element_blank(),
    legend.position = "none"
    )

fig.SL.45_60
```
# Grwoth in weight (Day45 and 60)
## 1. Random effect structures 
```{r}
lmer.W45a = lmer(W ~ Treatment * Day + (1|Tank) + (1|Clutch), data = MHWday45_60)
lmer.W45b = lmer(W ~ Treatment * Day + (1|Clutch/Tank), data = MHWday45_60)

AIC(lmer.W45a, lmer.W45b)
```

Choosing b. 

## 2. Covaraite exploration
```{r}
lm = lm(W ~ FemaleSL, data=MHWday45_60)
summary(lm)

lm = lm(W ~ FemaleFK, data=MHWday45_60)
summary(lm)

lm = lm(W ~ MaleSL, data=MHWday45_60)
summary(lm)

lm = lm(W ~ MaleFK, data=MHWday45_60)
summary(lm)
```

Correlation with Male SL. Testing model fit. 

```{r}
lmer.W45b = lmer(W ~ Treatment * Day + (1|Clutch/Tank), data = MHWday45_60)
lmer.W45c = lmer(W ~ Treatment * Day + MaleSL + (1|Clutch/Tank), data = MHWday45_60)

AIC(lmer.W45b, lmer.W45c)
```
Choosing b (without covariate).  

## 3. Final model 
```{r}
lmer.W45 = lmer(W ~ Treatment * Day + (1|Clutch/Tank), data = MHWday45_60)
```

## 4. Performance check
```{r}
performance::check_model(lmer.W45, check="homogeneity") 
performance::check_model(lmer.W45, check="outliers") 
performance::check_model(lmer.W45, check="qq", detrend = FALSE) 
performance::check_model(lmer.W45, check="normality") 
performance::check_model(lmer.W45, check="linearity") 
performance::check_model(lmer.W45, check="pp_check")
hist(residuals(lmer.W45), col="darkgray")
shapiro.test(residuals(lmer.W45))
qqline(resid(lmer.W45))
outlierTest(lmer.W45)
```
Outlier appears within the range of clutches. Kept. 

## 5. Outputs
### Summary statistic
```{r}
summary(lmer.W45)
```

### anova
```{r}
Anova(lmer.W45, type = 3)
```
### Emmeans
```{r}
EmmW.45_60 = (emmeans(lmer.W45, ~ Treatment * Day) %>% as.data.frame) 

EmmW.45_60
```

### Pairwise
#### Day 
```{r}
emm_day <- emmeans(lmer.W45, ~ Day)

pairs_day <- pairs(emm_day, adjust = "tukey") |> summary(infer = TRUE)
pairs_day
```


## 6. Figure
```{r}
fig.W.45_60 <- ggplot(EmmW.45_60, 
                        aes(x = Day, y = emmean, fill = Treatment)) + 
  geom_errorbar(
    aes(ymin = lower.CL, ymax = upper.CL),
    position = position_dodge(width = 0.6),     
    width = 0.1, colour = "black", linewidth = 0.3
  ) + 
  geom_point(
  aes(size = interaction(Treatment, Day)),  
  shape = 22, colour = "black", stroke = 0.3,
  position = position_dodge(width = 0.6)
) +
scale_size_manual(
  values = c(
    "Control.45" = 6,
    "Early exposure.45"  = 6,
    "Late exposure.45"  = 6,
    "Control.60" = 6,
    "Early exposure.60"  = 6,
    "Late exposure.60"  = 6
  ),
  guide = "none"  
) +
  theme_classic() +
  labs(
    y = "Weight (g)" 
  ) +
  scale_fill_manual(values = c("Control" = "#0072B2", 
                               "Early exposure" = "#F2C200",
                               "Late exposure" = "#D55E00")) +  
  scale_y_continuous(
  limits = c(0.1, 0.5),
  breaks = c(0.1, 0.3, 0.5),
  labels = c(0.1, 0.3, 0.5)) +
  theme(
    text = element_text(family = "Helvetica"),
    axis.title.x = element_blank(),
    axis.title = element_text(size = 14, face = "bold", colour = "black"),
    axis.text = element_text(size = 12, colour = "black"),
    axis.text.x  = element_blank(),
    legend.position = "none"
    )

fig.W.45_60
```

# Growth in Condition (Day 45 and 60)
## 1. Random effect structure
Exploring whether SL and W data should be logged. 
```{r}
ggplot(MHWday45_60, aes(x = SL, y = W)) + 
  geom_point() + 
  geom_smooth(method = "lm", col = "red") +
  ggtitle("Scatter Plot with Best Fit Line") +
  xlab("SL") +
  ylab("W")

ggplot(MHWday45_60, aes(x = SL, y = W)) +
  geom_point() +
  geom_smooth(method = "lm", formula = y ~ poly(x, 3), col = "blue") +
  ggtitle("Scatter Plot with Cubic Fit Line") +
   xlab("SL") +
  ylab("W")

ggplot(MHWday45_60, aes(x = LogSL, y = LogW)) + 
  geom_point() + 
  geom_smooth(method = "lm", col = "red") +
  ggtitle("Scatter Plot with Best Fit Line") +
  xlab("LogSL") +
  ylab("LogW")

ggplot(MHWday45_60, aes(x = LogSL, y = LogW)) +
  geom_point() +
  geom_smooth(method = "lm", formula = y ~ poly(x, 3), col = "blue") +
  ggtitle("Scatter Plot with Cubic Fit Line") +
   xlab("LogSL") +
  ylab("LogW")
```

Logging data to meet linearity. 

```{r}
MHWday45_60 <- MHWday45_60 %>%
  group_by(Day) %>%
  mutate(LogSL.centered = LogSL - mean(LogSL, na.rm = TRUE))
```

```{r}
lmer.Cond45a = lmer(LogW ~ Treatment * Day + LogSL.centered + (1|Tank) + (1|Clutch), data = MHWday45_60)
lmer.Cond45b = lmer(LogW ~ Treatment * Day + LogSL.centered + (1|Clutch/Tank), data = MHWday45_60)

AIC(lmer.Cond45a, lmer.Cond45b)
```

Choosing b. 

## 2. Covariate exploration 
Following weight, testing model fit including Male SL. 
```{r}
lmer.Cond45b = lmer(LogW ~ Treatment * Day + LogSL.centered + (1|Clutch/Tank), data = MHWday45_60)
lmer.Cond45c = lmer(LogW ~ Treatment * Day + LogSL.centered + LogMaleSL + (1|Clutch/Tank), data = MHWday45_60)


AIC(lmer.Cond45b, lmer.Cond45c)
```

Choosing b (without covariate). 

## 3. Final model 
```{r}
lmer.Cond45 = lmer(LogW ~ Treatment * Day + LogSL.centered + (1|Clutch/Tank), data = MHWday45_60)
```

## 4. Model performance
```{r}
performance::check_model(lmer.Cond45, check="homogeneity") 
performance::check_model(lmer.Cond45, check="outliers") 
performance::check_model(lmer.Cond45, check="qq", detrend = FALSE) 
performance::check_model(lmer.Cond45, check="normality") 
performance::check_model(lmer.Cond45, check="linearity") 
performance::check_model(lmer.Cond45, check="pp_check")
hist(residuals(lmer.Cond45), col="darkgray")
shapiro.test(residuals(lmer.Cond45))
qqline(resid(lmer.Cond45))
outlierTest(lmer.Cond45)
```
Both outlier appear to be within the expected range of a clutch. Kept. 

## 5. Outputs
### Summary statistics
```{r}
summary(lmer.Cond45)
```

### anova
```{r}
Anova(lmer.Cond45, type = 3)
```

### Emmeans
```{r}
EmmCond.45_60.Log = (emmeans(lmer.Cond45, ~ Treatment * Day) %>% as.data.frame) 

EmmCond.45_60 <- EmmCond.45_60.Log %>%
  mutate(
    emmean = exp(emmean),
    lower.CL = exp(lower.CL),
    upper.CL = exp(upper.CL)
  )

EmmCond.45_60
```

### Pairwise
#### Day 
```{r}
emm_day <- emmeans(lmer.Cond45, ~ Day)

pairs_day <- pairs(emm_day, adjust = "tukey") |> summary(infer = TRUE)
pairs_day
```

## 6. Figure
```{r}
fig.Cond.45_60 <- ggplot(EmmCond.45_60, 
                        aes(x = Day, y = emmean, fill = Treatment)) + 
  geom_errorbar(
    aes(ymin = lower.CL, ymax = upper.CL),
    position = position_dodge(width = 0.6),     
    width = 0.1, colour = "black", linewidth = 0.3
  ) + 
  geom_point(
  aes(size = interaction(Treatment, Day)),  
  shape = 22, colour = "black", stroke = 0.3,
  position = position_dodge(width = 0.6)
) +
scale_size_manual(
  values = c(
    "Control.45" = 6,
    "Early exposure.45"  = 6,
    "Late exposure.45"  = 6,
    "Control.60" = 6,
    "Early exposure.60"  = 6,
    "Late exposure.60"  = 6
  ),
  guide = "none"  
) +
  theme_classic() +
  labs(
    y = "Condition (g)" 
  ) +
  scale_fill_manual(values = c("Control" = "#0072B2", 
                               "Early exposure" = "#F2C200",
                               "Late exposure" = "#D55E00")) +   
  scale_y_continuous(
  limits = c(0.175, 0.4),
  breaks = c(0.2, 0.3, 0.4),
  labels = c(0.2, 0.3, 0.4)) +
  theme(
    text = element_text(family = "Helvetica"),
    axis.title.x = element_blank(),
    axis.title = element_text(size = 14, face = "bold", colour = "black"),
    axis.text = element_text(size = 12, colour = "black"),
    legend.position = "none"
    )

fig.Cond.45_60
```