R Markdown

library(lme4)
library(DHARMa)
library(ggplot2)
library(glmmTMB)
library(nlme)
library(MASS)
library(MuMIn)
library(vegan)
library(emmeans)
library(dplyr)
library(lmerTest)
library(car)
library(sjPlot)
library(patchwork)
library(ggpubr)
library(gridExtra)
library(readxl)
getwd()
## [1] "/Users/jc492396/Desktop/Publication/Honours_Ch2/MEPS/Data and Stats/4th submission"

Data import

F1 morphometric traits

Adults <- read_excel("ThermalHistoryCoralReefFish.xlsx", 
                     sheet = "Adult traits", 
                     trim_ws = TRUE) %>% 
          setNames(make.names(names(.), unique = TRUE))

Adults$Tank.ID = factor(Adults$Tank.ID)
Adults$Treatment = factor(Adults$Treatment)
Adults$Dev.treatment = factor(Adults$Dev.treatment)
Adults$Repro.treatment = factor(Adults$Repro.treatment)
Adults$Sex = factor(Adults$Sex)
Adults$Family = factor(Adults$Family)
Adults$Dev.tank = factor(Adults$Dev.tank)

str(Adults)
## tibble [86 × 15] (S3: tbl_df/tbl/data.frame)
##  $ Tank.ID        : Factor w/ 43 levels "1","6","9","11",..: 1 1 2 2 3 3 4 4 5 5 ...
##  $ Sump           : num [1:86] 1 1 1 1 1 1 1 1 6 6 ...
##  $ Fish.ID        : chr [1:86] "DHHj" "BHH3" "ACH" "BCHb" ...
##  $ Treatment      : Factor w/ 4 levels "Control-Control",..: 4 4 2 2 4 4 4 4 1 1 ...
##  $ Dev.treatment  : Factor w/ 2 levels "Control","Warm": 2 2 1 1 2 2 2 2 1 1 ...
##  $ Repro.treatment: Factor w/ 2 levels "Control","Warm": 2 2 2 2 2 2 2 2 1 1 ...
##  $ Sex            : Factor w/ 2 levels "F","M": 1 2 1 2 1 2 1 2 1 2 ...
##  $ Family         : Factor w/ 6 levels "A","B","C","D",..: 4 2 1 2 1 3 1 4 3 4 ...
##  $ Family.cross   : chr [1:86] "DB" "DB" "AB" "AB" ...
##  $ Dev.tank       : Factor w/ 54 levels "AC1","AC10","AC2",..: 51 22 2 49 14 33 9 44 27 34 ...
##  $ SL.cm          : num [1:86] 9.02 9.09 10.18 9.3 9.07 ...
##  $ Weight.g       : num [1:86] 32.2 32.2 42 33.1 33.2 ...
##  $ Fulton.K       : num [1:86] 4.39 4.28 3.98 4.11 4.44 ...
##  $ log10.SL       : num [1:86] 0.955 0.959 1.008 0.968 0.958 ...
##  $ log10.W        : num [1:86] 1.51 1.51 1.62 1.52 1.52 ...

Proportion of breeding

Breeding <- read_excel("ThermalHistoryCoralReefFish.xlsx", 
                     sheet = "Breeding proportions", 
                     trim_ws = TRUE) %>% 
          setNames(make.names(names(.), unique = TRUE))

Breeding$Development = factor(Breeding$Development)
Breeding$Post.maturity = factor(Breeding$Post.maturity)
Breeding$Treatment.combo = factor(Breeding$Treatment.combo)

str(Breeding)
## tibble [43 × 4] (S3: tbl_df/tbl/data.frame)
##  $ Development    : Factor w/ 2 levels "Control","Warm": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Post.maturity  : Factor w/ 2 levels "Control","Warm": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Treatment.combo: Factor w/ 4 levels "CC","CH","HC",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Reproduced     : chr [1:43] "N" "N" "N" "N" ...

Individual traits models

For each F1 adult traits, the information is presented as follows:
1. Final model and model performance:
2. Summary statistics:
Summary statistics are cited in Table S5 for each trait.
3. anova results:
Results from anova are cited in Table S2 for each trait.
4. Estimated marginal means (emmeans):
Emmeans and 95% confidence intervals are provided for each trait. These values are used to generate the corresponding plots.
5. Trait plots:
Plots for each traits are combined at the end for Figure 4.

Standard length

1. Final model and model performance

Model formula

lmer.SL = lmer(SL.cm  ~ Dev.treatment * Repro.treatment + (1|Family), data = Adults)

Model performance

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")

hist(residuals(lmer.SL), col="darkgray")

shapiro.test(residuals(lmer.SL))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(lmer.SL)
## W = 0.98279, p-value = 0.3094
qqnorm(resid(lmer.SL))
qqline(resid(lmer.SL))

outlierTest(lmer.SL)
## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
##    rstudent unadjusted p-value Bonferroni p
## 37 3.395807          0.0010733     0.092305

2. Summary statistics … cited in Table S5

summary(lmer.SL)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: SL.cm ~ Dev.treatment * Repro.treatment + (1 | Family)
##    Data: Adults
## 
## REML criterion at convergence: 141
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1568 -0.5768 -0.0844  0.5861  3.2219 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Family   (Intercept) 0.06677  0.2584  
##  Residual             0.25762  0.5076  
## Number of obs: 86, groups:  Family, 6
## 
## Fixed effects:
##                                       Estimate Std. Error       df t value
## (Intercept)                            8.97057    0.14600  9.34994  61.444
## Dev.treatmentWarm                     -0.07164    0.15908 79.30331  -0.450
## Repro.treatmentWarm                    0.17044    0.14259 78.29196   1.195
## Dev.treatmentWarm:Repro.treatmentWarm -0.44627    0.22737 78.37781  -1.963
##                                       Pr(>|t|)    
## (Intercept)                           1.65e-13 ***
## Dev.treatmentWarm                       0.6537    
## Repro.treatmentWarm                     0.2356    
## Dev.treatmentWarm:Repro.treatmentWarm   0.0532 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Dv.trW Rpr.tW
## Dv.trtmntWr -0.429              
## Rpr.trtmntW -0.478  0.462       
## Dv.trtW:R.W  0.300 -0.702 -0.636

3. anova … cited in Table S2

anova(lmer.SL)
## Type III Analysis of Variance Table with Satterthwaite's method
##                                Sum Sq Mean Sq NumDF  DenDF F value  Pr(>F)  
## Dev.treatment                 1.74106 1.74106     1 77.956  6.7583 0.01116 *
## Repro.treatment               0.05661 0.05661     1 76.884  0.2198 0.64055  
## Dev.treatment:Repro.treatment 0.99244 0.99244     1 78.378  3.8524 0.05322 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

4. emmeans and pairwise comparison

EmmSL = (emmeans(lmer.SL, ~ Dev.treatment * Repro.treatment) %>% as.data.frame) 
EmmSL
##  Dev.treatment Repro.treatment   emmean        SE    df lower.CL upper.CL
##  Control       Control         8.970572 0.1462098 10.40 8.646476 9.294667
##  Warm          Control         8.898932 0.1646588 14.95 8.547876 9.249989
##  Control       Warm            9.141010 0.1481658 10.51 8.813024 9.468995
##  Warm          Warm            8.623097 0.1670127 16.70 8.270244 8.975949
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95

Pairwise comparison is provided here as anova suggests for the interaction between developmental and post-maturation treatments.

emmeans(lmer.SL, pairwise ~ Dev.treatment * Repro.treatment, adjust = "tukey")$contrasts
##  contrast                       estimate    SE   df t.ratio p.value
##  Control Control - Warm Control   0.0716 0.160 79.6   0.447  0.9700
##  Control Control - Control Warm  -0.1704 0.143 78.7  -1.190  0.6347
##  Control Control - Warm Warm      0.3475 0.163 77.9   2.137  0.1506
##  Warm Control - Control Warm     -0.2421 0.157 78.0  -1.538  0.4202
##  Warm Control - Warm Warm         0.2758 0.176 77.8   1.569  0.4021
##  Control Warm - Warm Warm         0.5179 0.162 77.5   3.192  0.0108
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 4 estimates

5. Plot … Figure 2a

For all trait plots, post-maturation treatments are given in the x-axis, and developmental treatments shown in colors (blue = Control, orange = Warm).
Sample sizes for each of the 4 treatment combination are displayed to the right of each points, with raw values for each fish represented with jitterred points.
The x-axis (i.e., post-maturation treatment) is left unlabeled in individual plots to allow a common label to be later added when combining plots.

sample_sizes <- c(26, 16, 26, 18)

fig.2a <- ggplot(EmmSL, aes(x = Repro.treatment, y = emmean, colour = Dev.treatment)) + 
  geom_jitter(data = Adults, 
              aes(x = Repro.treatment, y = SL.cm, colour = Dev.treatment, 
                  group = interaction(Dev.treatment, Repro.treatment)),
              position = position_jitterdodge(jitter.width = 0.2, jitter.height = 0, dodge.width = 0.6), 
              size = 1.5, alpha = 0.2) + 
  geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL, group = interaction(Dev.treatment, Repro.treatment)), 
                position = position_dodge(0.6), linewidth = 0.4, colour = "black", width = 0.05) +
  geom_point(aes(group = interaction(Dev.treatment, Repro.treatment), fill = Dev.treatment),
             shape = 22, 
             position = position_dodge(0.6), size = 3, colour = "black", stroke = 0.5) +
  geom_text(aes(label = sample_sizes, group = interaction(Dev.treatment, Repro.treatment)), 
            position = position_dodge(0.6), vjust = -1, hjust = -0.75, size = 4, colour = "black", family = "Helvetica") +  
  theme_classic() +
  labs(x = NULL,
       y = "Standard length (cm)") + 
  scale_colour_manual(
    name = "Developmental treatment",  
    values = c("Control" = "blue", "Warm" = "orange")) + 
  scale_fill_manual(
    name = "Developmental treatment",  
    values = c("Control" = "blue", "Warm" = "orange")) + 
  guides(
    fill = guide_legend(
      title.position = "top",  
      title.hjust = 0.5,       
      nrow = 2,                
      direction = "horizontal" 
    )
  ) +
  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",                 
    legend.title = element_text(size = 12, face = "bold", colour = "black"),  
    legend.text = element_text(size = 12, colour = "black"),
    legend.background = element_rect(fill = "white", colour = "black"),  
    legend.box.background = element_rect(colour = "black", linewidth = 0.1)
  )

fig.2a

Weight

1. Final model and model performance

Model formula

lmer.W = lmer(Weight.g ~ Dev.treatment * Repro.treatment + (1|Family) , data = Adults)

Model performance

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))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(lmer.W)
## W = 0.93459, p-value = 0.0002964
qqnorm(resid(lmer.W))
qqline(resid(lmer.W))

outlierTest(lmer.W)
##    rstudent unadjusted p-value Bonferroni p
## 37 4.369868         3.7434e-05    0.0032193

Fish 37 was just a large fish. Not removed as outlier.

2. Summary statistics … cited in Table S5

summary(lmer.W)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Weight.g ~ Dev.treatment * Repro.treatment + (1 | Family)
##    Data: Adults
## 
## REML criterion at convergence: 532.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6616 -0.6754 -0.1629  0.6317  4.1460 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Family   (Intercept)  7.925   2.815   
##  Residual             30.468   5.520   
## Number of obs: 86, groups:  Family, 6
## 
## Fixed effects:
##                                       Estimate Std. Error      df t value
## (Intercept)                            31.6066     1.5892  8.9722  19.888
## Dev.treatmentWarm                      -0.9827     1.7301 79.1797  -0.568
## Repro.treatmentWarm                    -0.5700     1.5507 78.1277  -0.368
## Dev.treatmentWarm:Repro.treatmentWarm  -2.5631     2.4727 78.2167  -1.037
##                                       Pr(>|t|)    
## (Intercept)                           9.94e-09 ***
## Dev.treatmentWarm                        0.572    
## Repro.treatmentWarm                      0.714    
## Dev.treatmentWarm:Repro.treatmentWarm    0.303    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Dv.trW Rpr.tW
## Dv.trtmntWr -0.428              
## Rpr.trtmntW -0.478  0.462       
## Dv.trtW:R.W  0.300 -0.702 -0.636

3. anova … cited in Table S2

anova(lmer.W)
## Type III Analysis of Variance Table with Satterthwaite's method
##                                Sum Sq Mean Sq NumDF  DenDF F value  Pr(>F)  
## Dev.treatment                 102.725 102.725     1 77.778  3.3716 0.07015 .
## Repro.treatment                69.884  69.884     1 76.665  2.2937 0.13401  
## Dev.treatment:Repro.treatment  32.735  32.735     1 78.217  1.0744 0.30314  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

4. Emmeans

EmmW = (emmeans(lmer.W, ~ Dev.treatment * Repro.treatment) %>% as.data.frame) 
EmmW
##  Dev.treatment Repro.treatment   emmean       SE    df lower.CL upper.CL
##  Control       Control         31.60661 1.591552 10.38 28.07794 35.13529
##  Warm          Control         30.62387 1.792045 14.92 26.80254 34.44521
##  Control       Warm            31.03662 1.612819 10.49 27.46570 34.60754
##  Warm          Warm            27.49080 1.817606 16.66 23.65002 31.33159
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95

5. Plot … Figure 2b

fig.2b <- ggplot(EmmW, aes(x = Repro.treatment, y = emmean, colour = Dev.treatment)) + 
  geom_jitter(data = Adults, 
              aes(x = Repro.treatment, y = Weight.g, colour = Dev.treatment, 
                  group = interaction(Dev.treatment, Repro.treatment)),
              position = position_jitterdodge(jitter.width = 0.2, jitter.height = 0, dodge.width = 0.6), 
              size = 1.5, alpha = 0.2) + 
  geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL, group = interaction(Dev.treatment, Repro.treatment)), 
                position = position_dodge(0.6), linewidth = 0.4, colour = "black", width = 0.05) + 
  geom_point(aes(group = interaction(Dev.treatment, Repro.treatment), fill = Dev.treatment),
             shape = 22, 
             position = position_dodge(0.6), size = 3, colour = "black", stroke = 0.5) +
  theme_classic() +
  labs(x = NULL,
       y = "Weight (g)") + 
  scale_fill_manual(values = c("Control" = "blue", "Warm" = "orange")) +  
  scale_colour_manual(values = c("Control" = "blue", "Warm" = "orange")) + 
  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.2b

Condition (W for given SL)

1. Final model and model performance

Model formula

lmer.Cond = lmer(Weight.g ~ Dev.treatment * Repro.treatment + (1|Family) + SL.cm, data = Adults)

Performance check

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))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(lmer.Cond)
## W = 0.98113, p-value = 0.2418
qqnorm(resid(lmer.Cond))
qqline(resid(lmer.Cond))

outlierTest(lmer.Cond)
## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
##    rstudent unadjusted p-value Bonferroni p
## 37  3.14584          0.0023446      0.20163

2. Summary statistics … cited in Table S5

summary(lmer.Cond)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Weight.g ~ Dev.treatment * Repro.treatment + (1 | Family) + SL.cm
##    Data: Adults
## 
## REML criterion at convergence: 402.9
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.18392 -0.62982 -0.05598  0.52896  2.79294 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Family   (Intercept) 1.034    1.017   
##  Residual             6.548    2.559   
## Number of obs: 86, groups:  Family, 6
## 
## Fixed effects:
##                                       Estimate Std. Error       df t value
## (Intercept)                           -55.2663     4.9772  79.8469 -11.104
## Dev.treatmentWarm                      -0.2953     0.7999  79.1758  -0.369
## Repro.treatmentWarm                    -2.2156     0.7239  77.4951  -3.060
## SL.cm                                   9.6840     0.5499  80.3915  17.611
## Dev.treatmentWarm:Repro.treatmentWarm   1.7745     1.1705  77.8895   1.516
##                                       Pr(>|t|)    
## (Intercept)                            < 2e-16 ***
## Dev.treatmentWarm                      0.71298    
## Repro.treatmentWarm                    0.00304 ** 
## SL.cm                                  < 2e-16 ***
## Dev.treatmentWarm:Repro.treatmentWarm  0.13357    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Dv.trW Rpr.tW SL.cm 
## Dv.trtmntWr -0.111                     
## Rpr.trtmntW  0.064  0.449              
## SL.cm       -0.991  0.048 -0.136       
## Dv.trtW:R.W -0.168 -0.673 -0.643  0.214

3. anova … cited in Table S2

anova(lmer.Cond)
## Type III Analysis of Variance Table with Satterthwaite's method
##                                Sum Sq Mean Sq NumDF  DenDF  F value  Pr(>F)    
## Dev.treatment                    6.51    6.51     1 77.889   0.9949 0.32165    
## Repro.treatment                 35.92   35.92     1 76.697   5.4858 0.02177 *  
## SL.cm                         2030.93 2030.93     1 80.391 310.1399 < 2e-16 ***
## Dev.treatment:Repro.treatment   15.05   15.05     1 77.890   2.2983 0.13357    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

4. emmeans

EmmCond = (emmeans(lmer.Cond, ~ Dev.treatment * Repro.treatment) %>% as.data.frame)
EmmCond
##  Dev.treatment Repro.treatment   emmean        SE    df lower.CL upper.CL
##  Control       Control         31.65649 0.6580959 12.84 30.23300 33.07998
##  Warm          Control         31.36118 0.7599532 18.89 29.76997 32.95239
##  Control       Warm            29.44093 0.6732899 13.49 27.99170 30.89016
##  Warm          Warm            30.92017 0.7987959 23.62 29.27013 32.57022
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95

5. Plots … Figure 2c

The fitted values from the model is extracted first to plot the jitterred points.

Cond_fitted <- predict(lmer.Cond, newdata = Adults, re.form = NULL)
Adults <- cbind(Adults[, 1:13], Cond_fitted, Adults[, 14:ncol(Adults)])
colnames(Adults)[14] <- "Cond_fitted"
fig.2c <- ggplot(EmmCond, aes(x = Repro.treatment, y = emmean, colour = Dev.treatment)) + 
  geom_jitter(data = Adults, 
              aes(x = Repro.treatment, y = Cond_fitted, colour = Dev.treatment, 
                  group = interaction(Dev.treatment, Repro.treatment)),
              position = position_jitterdodge(jitter.width = 0.2, jitter.height = 0, dodge.width = 0.6), 
              size = 1.5, alpha = 0.2) +  
  geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL, group = interaction(Dev.treatment, Repro.treatment)), 
                position = position_dodge(0.6), linewidth = 0.4, colour = "black", width = 0.05) +  
  geom_point(aes(group = interaction(Dev.treatment, Repro.treatment), fill = Dev.treatment),
             shape = 22, 
             position = position_dodge(0.6), size = 2.8, colour = "black", stroke = 0.5) +
  theme_classic() +
  labs(x = NULL,
       y = "Condition") + 
  scale_fill_manual(values = c("Control" = "blue", "Warm" = "orange")) +  
  scale_colour_manual(values = c("Control" = "blue", "Warm" = "orange")) + 
  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.2c

Fulton’s K

1. Final model and model performance

Model formula

lmer.FK = lmer(Fulton.K ~ Dev.treatment * Repro.treatment + (1|Family) , data = Adults)

Performance checks

performance::check_model(lmer.FK, check="homogeneity") 

performance::check_model(lmer.FK, check="outliers")

performance::check_model(lmer.FK, check="qq", detrend = FALSE) 

performance::check_model(lmer.FK, check="normality") 

performance::check_model(lmer.FK, check="linearity") 

performance::check_model(lmer.FK, check="pp_check")

hist(residuals(lmer.FK), col="darkgray")

shapiro.test(residuals(lmer.FK))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(lmer.FK)
## W = 0.9714, p-value = 0.05319
qqnorm(resid(lmer.FK))
qqline(resid(lmer.FK))

outlierTest(lmer.FK)
## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
##    rstudent unadjusted p-value Bonferroni p
## 73 3.447477          0.0009098     0.078243

2. Summary statistics … cited in Table S5

summary(lmer.FK)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Fulton.K ~ Dev.treatment * Repro.treatment + (1 | Family)
##    Data: Adults
## 
## REML criterion at convergence: 81.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.2340 -0.6127 -0.1535  0.4967  3.3017 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Family   (Intercept) 0.01544  0.1243  
##  Residual             0.12879  0.3589  
## Number of obs: 86, groups:  Family, 6
## 
## Fixed effects:
##                                       Estimate Std. Error       df t value
## (Intercept)                            4.34114    0.08741 14.55042  49.662
## Dev.treatmentWarm                     -0.06280    0.11180 80.35515  -0.562
## Repro.treatmentWarm                   -0.30265    0.10045 79.15365  -3.013
## Dev.treatmentWarm:Repro.treatmentWarm  0.26862    0.16014 79.31948   1.677
##                                       Pr(>|t|)    
## (Intercept)                            < 2e-16 ***
## Dev.treatmentWarm                      0.57585    
## Repro.treatmentWarm                    0.00347 ** 
## Dev.treatmentWarm:Repro.treatmentWarm  0.09740 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Dv.trW Rpr.tW
## Dv.trtmntWr -0.509              
## Rpr.trtmntW -0.566  0.459       
## Dv.trtW:R.W  0.355 -0.699 -0.634

3. anova … cited in Table S2

anova(lmer.FK)
## Type III Analysis of Variance Table with Satterthwaite's method
##                                Sum Sq Mean Sq NumDF  DenDF F value  Pr(>F)  
## Dev.treatment                 0.10308 0.10308     1 78.880  0.8004 0.37371  
## Repro.treatment               0.57866 0.57866     1 77.487  4.4931 0.03723 *
## Dev.treatment:Repro.treatment 0.36237 0.36237     1 79.319  2.8137 0.09740 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

4. Emmeans

EmmFK = (emmeans(lmer.FK, ~ Dev.treatment * Repro.treatment) %>% as.data.frame) 
EmmFK
##  Dev.treatment Repro.treatment   emmean         SE    df lower.CL upper.CL
##  Control       Control         4.341144 0.08774688 15.14 4.154262 4.528026
##  Warm          Control         4.278341 0.10222516 22.14 4.066415 4.490267
##  Control       Warm            4.038490 0.08906114 14.80 3.848432 4.228547
##  Warm          Warm            4.244308 0.10449092 26.59 4.029754 4.458862
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95

5. Plot … Figure S1

sample_sizes <- c(26, 16, 26, 18)

fig.S1 <- ggplot(EmmFK, aes(x = Repro.treatment, y = emmean, colour = Dev.treatment)) + 
  geom_jitter(data = Adults, 
              aes(x = Repro.treatment, y = Fulton.K, colour = Dev.treatment, 
                  group = interaction(Dev.treatment, Repro.treatment)),
              position = position_jitterdodge(jitter.width = 0.2, jitter.height = 0, dodge.width = 0.6), 
              size = 1.5, alpha = 0.2) + 
  geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL, group = interaction(Dev.treatment, Repro.treatment)), 
                position = position_dodge(0.6), linewidth = 0.4, colour = "black", width = 0.05) +  
  geom_point(aes(group = interaction(Dev.treatment, Repro.treatment), fill = Dev.treatment),
             shape = 22, 
             position = position_dodge(0.6), size = 3, colour = "black", stroke = 0.5) +
  geom_text(aes(label = sample_sizes, group = interaction(Dev.treatment, Repro.treatment)), 
            position = position_dodge(0.6), vjust = -1, hjust = -0.3, size = 4, colour = "black", family = "Helvetica") +  
  theme_classic() +
  labs(x = "Post-maturation treatment",
     y = "Fulton's K" ) +
  scale_colour_manual(
    name = "Developmental treatment",  
    values = c("Control" = "blue", "Warm" = "orange")) + 
  scale_fill_manual(
    name = "Developmental treatment",  
    values = c("Control" = "blue", "Warm" = "orange")) + 
  guides(
    fill = guide_legend(
      title = "Development", 
      title.position = "top",    
      label.hjust = 0.5,             
      direction = "horizontal",     
      nrow = 2                      
    ),
    colour = "none"
  ) +
  theme(
    text = element_text(family = "Helvetica"),  
    axis.title.x = element_text(size = 14, face = "bold", colour = "black"),
    axis.title.y = element_text(size = 14, face = "bold", colour = "black"),
    axis.text = element_text(size = 12, colour = "black"),
    legend.position = c(0.775, 0.85),                 
    legend.title = element_text(size = 12, face = "bold", colour = "black"),  
    legend.text = element_text(size = 12, colour = "black"),
    legend.background = element_rect(fill = "white", colour = "black"),  
    legend.box.background = element_rect(colour = "black", linewidth = 0.1)
  )

fig.S1

ggsave(
  filename = "/Users/jc492396/Desktop/Publication/Honours_Ch2/MEPS/Data and Stats/4th submission/Figures/Figure S1.pdf", 
  plot = fig.S1, 
  width = 89, 
  height = 100, 
  units = "mm", 
  dpi = 600
)

Proporation of breeding

The following analyses assess treatment effects on breeding proportion.
The X-squared values are cited in Results to discuss the influence of developmental and/or post-maturation treatment on breeding occurrence.

table(Breeding$Treatment.combo, Breeding$Reproduced)
##     
##      N Y
##   CC 7 6
##   CH 7 6
##   HC 4 5
##   HH 5 3

Does treatment (both during development and/or post-maturation) influence the occurrence of breeding?

chisq.test(table(Breeding$Treatment.combo, Breeding$Reproduced))
## 
##  Pearson's Chi-squared test
## 
## data:  table(Breeding$Treatment.combo, Breeding$Reproduced)
## X-squared = 0.55837, df = 3, p-value = 0.9059

Does warming during development or post-maturation influence the occurrence of breeding?

chisq.test(table(Breeding$Development, Breeding$Reproduced))
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  table(Breeding$Development, Breeding$Reproduced)
## X-squared = 0, df = 1, p-value = 1
chisq.test(table(Breeding$Post.maturity, Breeding$Reproduced))
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  table(Breeding$Post.maturity, Breeding$Reproduced)
## X-squared = 0.026759, df = 1, p-value = 0.8701

Combining plots for Figure 2

fig.2_legend <- ggplotGrob(fig.2a + theme(legend.position = "top"))
legend <- fig.2_legend$grobs[[which(sapply(fig.2_legend$grobs, function(x) x$name) == "guide-box")]]

fig.Adult <- ggarrange(
  fig.2a, NULL, fig.2b, fig.2c,
  labels = c("a", "", "b", "c"),     
  ncol = 2, nrow = 2,                 
  label.x = 0.01, label.y = 1.0,     
  font.label = list(size = 14, face = "bold"),  
  align = "v",                         
  heights = c(1, 1)                    
)

fig.Adult_with_legend <- ggarrange(
  legend, fig.Adult,
  ncol = 1, heights = c(0.1, 0.9)  
)

fig.Adult_with_legend <- annotate_figure(
  fig.Adult_with_legend,
  bottom = text_grob("Post-maturation treatment", 
                     size = 16, face = "bold", family = "Helvetica")
)

ggsave(
  filename = "/Users/jc492396/Desktop/Publication/Honours_Ch2/MEPS/Data and Stats/4th submission/Figures/Figure 2.pdf", 
  plot = fig.Adult_with_legend, 
  width = 169, 
  height = 180, 
  units = "mm", 
  dpi = 300
)