Thermal history influences reproduction and offspring quality in a coral reef fish: Newly hatched traits
Y Yasutake
March 2022, updated January 2025
library(lme4)
library(DHARMa)
library(ggplot2)
library(glmmTMB)
library(nlme)
library(MASS)
library(MuMIn)
library(vegan)
library(emmeans)
library(lmerTest)
library(car)
library(sjPlot)
library(patchwork)
library(ggpubr)
library(gridExtra)
library(dplyr)
library(plyr)
library(readxl)getwd()## [1] "/Users/jc492396/Desktop/Publication/Honours_Ch2/MEPS/Data and Stats/4th submission"Data loading
Juvenile traits at hatching
Juv <- read_excel("ThermalHistoryCoralReefFish.xlsx",
sheet = "Juveniles at hatching",
trim_ws = TRUE) %>%
setNames(make.names(names(.), unique = TRUE))
Juv$Sump = factor (Juv$Sump)
Juv$Treatment = factor (Juv$Treatment)
Juv$Dev.treatment = factor (Juv$Dev.treatment)
Juv$Repro.treatment = factor (Juv$Repro.treatment)
Juv$Clutch.replicate = factor (Juv$Clutch.replicate)
Juv$Fish.replicate = factor (Juv$Fish.replicate)
Juv$Male.family = factor (Juv$Male.family)
Juv$Female.family = factor (Juv$Female.family)
str(Juv)## tibble [387 × 27] (S3: tbl_df/tbl/data.frame)
## $ Tank.ID : num [1:387] 9 9 9 9 9 9 9 9 9 9 ...
## $ Sump : Factor w/ 6 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Treatment : Factor w/ 4 levels "Control-Control",..: 4 4 4 4 4 4 4 4 4 4 ...
## $ Dev.treatment : Factor w/ 2 levels "Control","Warm": 2 2 2 2 2 2 2 2 2 2 ...
## $ Repro.treatment : Factor w/ 2 levels "Control","Warm": 2 2 2 2 2 2 2 2 2 2 ...
## $ Clutch.replicate : Factor w/ 20 levels "Control_Control_1",..: 18 18 18 18 18 18 18 18 18 18 ...
## $ Fish.replicate : Factor w/ 387 levels "Control_Control_1_1",..: 327 338 340 341 342 343 344 345 346 328 ...
## $ Female.ID : chr [1:387] "AHH" "AHH" "AHH" "AHH" ...
## $ Male.ID : chr [1:387] "CHH9" "CHH9" "CHH9" "CHH9" ...
## $ Female.family : Factor w/ 5 levels "A","B","C","D",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Male.family : Factor w/ 5 levels "A","C","D","E",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ Family.combo : chr [1:387] "AC" "AC" "AC" "AC" ...
## $ Male.SL..cm. : num [1:387] 9.48 9.48 9.48 9.48 9.48 ...
## $ Male.weight..g. : num [1:387] 31.9 31.9 31.9 31.9 31.9 31.9 31.9 31.9 31.9 31.9 ...
## $ Male.condition : num [1:387] 3.75 3.75 3.75 3.75 3.75 ...
## $ Female.SL..cm. : num [1:387] 9.07 9.07 9.07 9.07 9.07 ...
## $ Female.weight..g. : num [1:387] 33.2 33.2 33.2 33.2 33.2 33.2 33.2 33.2 33.2 33.2 ...
## $ Female.condition : num [1:387] 4.44 4.44 4.44 4.44 4.44 ...
## $ Avg.egg.area : num [1:387] 5.03 5.03 5.03 5.03 5.03 ...
## $ Juvenile.weight..mg. : num [1:387] 3.3 3.3 3.3 3.5 3.5 3.6 3.6 3.6 3.6 3.7 ...
## $ Juvenile.SL : num [1:387] 4.88 4.89 4.93 5.06 5.1 ...
## $ Juvenile.yolk.area..mm.2.: num [1:387] 1.45 1.5 1.71 1.84 1.78 ...
## $ Fulton.K : num [1:387] 2.84 2.81 2.75 2.7 2.64 ...
## $ log.W : num [1:387] 0.519 0.519 0.519 0.544 0.544 ...
## $ log.SL : num [1:387] 0.688 0.69 0.693 0.704 0.707 ...
## $ log.FK : num [1:387] 0.453 0.449 0.439 0.432 0.422 ...
## $ Notes : chr [1:387] NA NA NA NA ...Correlation between maternal traits
Parental <- read_excel("ThermalHistoryCoralReefFish.xlsx",
sheet = "Correlation to parent traits",
trim_ws = TRUE) %>%
setNames(make.names(names(.), unique = TRUE))
Parental$Sump = factor(Parental$Sump)
Parental$Tank.ID = factor(Parental$Tank.ID)
Parental$Treatment = factor(Parental$Treatment)
Parental$Dev.treatment = factor(Parental$Dev.treatment)
Parental$Repro.treatment = factor(Parental$Repro.treatment)
Parental$Pair.replicate = factor(Parental$Pair.replicate)
Parental$Female.family = factor(Parental$Male.family)
Parental$Male.family = factor(Parental$Male.family)
Parental$Male.family = factor(Parental$Male.family)
str(Parental)## tibble [20 × 25] (S3: tbl_df/tbl/data.frame)
## $ Tank.ID : Factor w/ 20 levels "9","25","42",..: 11 14 15 16 18 20 2 5 6 7 ...
## $ Sump : Factor w/ 6 levels "1","2","3","4",..: 4 6 5 5 5 5 1 3 3 3 ...
## $ Treatment : Factor w/ 4 levels "Control-Control",..: 1 1 1 1 1 1 2 2 2 2 ...
## $ Dev.treatment : Factor w/ 2 levels "Control","Warm": 1 1 1 1 1 1 1 1 1 1 ...
## $ Repro.treatment : Factor w/ 2 levels "Control","Warm": 1 1 1 1 1 1 2 2 2 2 ...
## $ Pair.replicate : Factor w/ 20 levels "Control-Control-1",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ Female.ID : chr [1:20] "ACC4" "ACC3" "ECC" "ECCb" ...
## $ Female.family : Factor w/ 5 levels "A","C","D","E",..: 5 3 1 5 2 1 4 3 4 3 ...
## $ SL...Female..cm. : num [1:20] 9.14 8.91 9.27 7.8 8.78 ...
## $ Weight...Female..g. : num [1:20] 31.3 35.9 33 26.6 32.6 ...
## $ Condition...Female : num [1:20] 4.09 5.07 4.14 5.61 4.83 ...
## $ Male.ID : chr [1:20] "FCC" "DCC" "ACC2" "FCC" ...
## $ Male.family : Factor w/ 5 levels "A","C","D","E",..: 5 3 1 5 2 1 4 3 4 3 ...
## $ Family.combo : chr [1:20] "AC" "AC" "AC" "AC" ...
## $ SL...Male..cm. : num [1:20] 8.65 9.65 9.45 7.74 8.6 ...
## $ Weight...Male..g. : num [1:20] 22.6 41.1 36.9 19.7 23.6 ...
## $ Condition...Male : num [1:20] 3.5 4.58 4.37 4.25 3.71 ...
## $ Total.number.of.eggs : num [1:20] 711 940 604 491 992 606 778 913 883 782 ...
## $ Number.of.clutches.per.pair : num [1:20] 3 3 2 3 3 2 3 3 3 3 ...
## $ Clutch.size..1st. : num [1:20] 198 324 309 123 349 302 216 323 364 278 ...
## $ Average.of.Egg.area..mm2. : num [1:20] 4.95 5.88 5.28 5.26 4.71 ...
## $ Average.of.Juvenile.weight..mg. : num [1:20] 3.85 4.03 3.55 3.88 4.11 ...
## $ Average.of.Juvenile.SL..mm. : num [1:20] 5.22 5.72 5.21 5.27 4.89 ...
## $ Average.of.Juvenile.yolk.area..mm.2.: num [1:20] 1.33 1.62 1.33 1.57 1.43 ...
## $ Average.of.Fulton.s.K : num [1:20] 2.71 2.15 2.51 2.65 3.52 ...Individual traits
For each newly hatched trait, the information is presented as
follows:
1. Maternal effects and covariate exploration:
A linear regression is performed between each newly hatched trait and
maternal physical variables.
A maternal variable is included in the final model if it improves model
fit based on AIC values.
2. Final model and model performance:
3. Summary statistics:
Summary statistics are cited in Table S7 for each trait.
4. anova results:
Results from anova are cited in Table S4 for each trait.
5. Estimated marginal means (emmeans):
Emmeans and 95% confidence intervals are provided for each trait. These
values are used to generate the corresponding plots.
6. Trait plots:
Plots for each traits are provided, which are combined at the end of the
markdown to create Figure 4.
Standard length
1. Maternal effects and covariate exploration
SLgroup_meanJSL = ddply(Juv, "Female.SL..cm.", summarise, grp.mean=mean(Juvenile.SL))
Wgroup_meanJSL = ddply(Juv, "Female.weight..g.", summarise, grp.mean=mean(Juvenile.SL))
Condgroup_meanJSL = ddply(Juv, "Female.condition", summarise, grp.mean=mean(Juvenile.SL))
Eggareagroup_meanJSL = ddply(Juv, "Avg.egg.area", summarise, grp.mean=mean(Juvenile.SL))Correlation between Female SL
lm.1 = lm(grp.mean ~ Female.SL..cm., data=SLgroup_meanJSL)
summary(lm.1)##
## Call:
## lm(formula = grp.mean ~ Female.SL..cm., data = SLgroup_meanJSL)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.43293 -0.11326 -0.00078 0.06775 0.50561
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.34847 1.10091 4.858 0.000126 ***
## Female.SL..cm. -0.01504 0.12373 -0.122 0.904606
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2505 on 18 degrees of freedom
## Multiple R-squared: 0.0008201, Adjusted R-squared: -0.05469
## F-statistic: 0.01477 on 1 and 18 DF, p-value: 0.9046Correaltion between Female W
lm.2 = lm(grp.mean ~ Female.weight..g., data=Wgroup_meanJSL)
summary(lm.2)##
## Call:
## lm(formula = grp.mean ~ Female.weight..g., data = Wgroup_meanJSL)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.49689 -0.12054 0.00414 0.08445 0.47312
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.063710 0.369795 13.693 5.86e-11 ***
## Female.weight..g. 0.005106 0.012353 0.413 0.684
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2494 on 18 degrees of freedom
## Multiple R-squared: 0.009405, Adjusted R-squared: -0.04563
## F-statistic: 0.1709 on 1 and 18 DF, p-value: 0.6842Correlation between Female FK
lm.3 = lm(grp.mean ~ Female.condition, data=Condgroup_meanJSL)
summary(lm.3)##
## Call:
## lm(formula = grp.mean ~ Female.condition, data = Condgroup_meanJSL)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.43467 -0.08638 -0.00717 0.08737 0.43031
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.80481 0.45732 10.506 4.15e-09 ***
## Female.condition 0.09732 0.10776 0.903 0.378
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2451 on 18 degrees of freedom
## Multiple R-squared: 0.04335, Adjusted R-squared: -0.009802
## F-statistic: 0.8156 on 1 and 18 DF, p-value: 0.3784No maternal traits for potential covariate.
2. Final model and model performance
Model formula
lmer.JSL = lmer(Juvenile.SL~ Dev.treatment * Repro.treatment + (1 |Sump/Female.ID), data = Juv)Performance checks
performance::check_model(lmer.JSL, check="homogeneity") 
performance::check_model(lmer.JSL, check="outliers")
performance::check_model(lmer.JSL, check="qq", detrend = FALSE) 
performance::check_model(lmer.JSL, check="normality") 
performance::check_model(lmer.JSL, check="linearity") 
performance::check_model(lmer.JSL, check="pp_check")
hist(residuals(lmer.JSL), col="darkgray")
outlierTest(lmer.JSL)## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
## rstudent unadjusted p-value Bonferroni p
## 198 -3.414145 0.00070895 0.274363. Summary statistics … cited in Table S7
summary(lmer.JSL)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Juvenile.SL ~ Dev.treatment * Repro.treatment + (1 | Sump/Female.ID)
## Data: Juv
##
## REML criterion at convergence: 12.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3304 -0.5618 0.0243 0.6800 2.7176
##
## Random effects:
## Groups Name Variance Std.Dev.
## Female.ID:Sump (Intercept) 0.061031 0.24705
## Sump (Intercept) 0.004319 0.06572
## Residual 0.050409 0.22452
## Number of obs: 387, groups: Female.ID:Sump, 20; Sump, 6
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 5.25169 0.11245 4.40934 46.703
## Dev.treatmentWarm 0.03070 0.15431 14.13894 0.199
## Repro.treatmentWarm -0.05194 0.15772 4.97389 -0.329
## Dev.treatmentWarm:Repro.treatmentWarm -0.14789 0.23709 14.30390 -0.624
## Pr(>|t|)
## (Intercept) 4.11e-07 ***
## Dev.treatmentWarm 0.845
## Repro.treatmentWarm 0.755
## Dev.treatmentWarm:Repro.treatmentWarm 0.543
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Dv.trW Rpr.tW
## Dv.trtmntWr -0.639
## Rpr.trtmntW -0.713 0.456
## Dv.trtW:R.W 0.416 -0.651 -0.5804. anova … cited in the Table S4
anova(lmer.JSL)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Dev.treatment 0.006708 0.006708 1 14.3039 0.1331 0.7206
## Repro.treatment 0.046344 0.046344 1 2.5837 0.9194 0.4185
## Dev.treatment:Repro.treatment 0.019615 0.019615 1 14.3039 0.3891 0.54265. emmeans
JSLEMM = (emmeans(lmer.JSL, ~ Dev.treatment*Repro.treatment) %>% as.data.frame)
JSLEMM## Dev.treatment Repro.treatment emmean SE df lower.CL upper.CL
## Control Control 5.251688 0.1311467 4.72 4.908377 5.595000
## Warm Control 5.282393 0.1219526 8.46 5.003829 5.560956
## Control Warm 5.199747 0.1166017 6.00 4.914483 5.485010
## Warm Warm 5.082558 0.1690999 10.52 4.708280 5.456835
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95emmeans(lmer.JSL, pairwise ~ Dev.treatment * Repro.treatment, adjust="tukey")## $emmeans
## Dev.treatment Repro.treatment emmean SE df lower.CL upper.CL
## Control Control 5.25 0.131 4.72 4.91 5.59
## Warm Control 5.28 0.122 8.46 5.00 5.56
## Control Warm 5.20 0.117 6.00 4.91 5.49
## Warm Warm 5.08 0.169 10.52 4.71 5.46
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## Control Control - Warm Control -0.0307 0.165 14.25 -0.186 0.9976
## Control Control - Control Warm 0.0519 0.175 5.30 0.296 0.9900
## Control Control - Warm Warm 0.1691 0.214 7.92 0.790 0.8569
## Warm Control - Control Warm 0.0826 0.169 7.21 0.490 0.9590
## Warm Control - Warm Warm 0.1998 0.208 9.70 0.958 0.7752
## Control Warm - Warm Warm 0.1172 0.189 14.50 0.621 0.9237
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 4 estimatesPlot … Figure 4a
For all trait plots, post-maturation treatments are given in the
x-axis, developmental treatments shown with 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 newly hatched 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_size_Juv <- c(94, 112, 61, 120)
fig.4a <- ggplot(JSLEMM, aes(x = Repro.treatment, y = emmean, colour = Dev.treatment)) +
geom_jitter(data = Juv,
aes(x = Repro.treatment, y = Juvenile.SL, 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 = 4, colour = "black", stroke = 0.5) +
geom_text(aes(label = sample_size_Juv, 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 (mm)") +
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 = "left",
label.hjust = 0,
direction = "horizontal",
nrow = 1
)
) +
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.4a
Weight
1. Maternal effects and covariate exploration
SLgroup_meanJW = ddply(Juv, "Female.SL..cm.", summarise, grp.mean=mean(Juvenile.weight..mg.))
Wgroup_meanJW = ddply(Juv, "Female.weight..g.", summarise, grp.mean=mean(Juvenile.weight..mg.))
Condgroup_meanJW = ddply(Juv, "Female.condition", summarise, grp.mean=mean(Juvenile.weight..mg.))
Egggroup_meanJW = ddply(Juv, "Avg.egg.area", summarise, grp.mean=mean(Juvenile.weight..mg.))Correlation with Female SL
lm.1 = lm(grp.mean ~ Female.SL..cm., data=SLgroup_meanJW)
summary(lm.1)##
## Call:
## lm(formula = grp.mean ~ Female.SL..cm., data = SLgroup_meanJW)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.49788 -0.27109 -0.01369 0.13215 0.59281
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.6310 1.6702 3.970 0.000989 ***
## Female.SL..cm. -0.2975 0.1889 -1.575 0.133631
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3313 on 17 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.1274, Adjusted R-squared: 0.07604
## F-statistic: 2.481 on 1 and 17 DF, p-value: 0.1336Correlation with Female W
lm.2 = lm(grp.mean ~ Female.weight..g., data=Wgroup_meanJW)
summary(lm.2)##
## Call:
## lm(formula = grp.mean ~ Female.weight..g., data = Wgroup_meanJW)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.6187 -0.2091 -0.0900 0.2132 0.6491
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.83977 0.55493 8.721 1.1e-07 ***
## Female.weight..g. -0.02877 0.01890 -1.523 0.146
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3327 on 17 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.12, Adjusted R-squared: 0.06826
## F-statistic: 2.319 on 1 and 17 DF, p-value: 0.1462Correlation with Female FK
lm.3 = lm(grp.mean ~ Female.condition, data=Condgroup_meanJW)
summary(lm.3)##
## Call:
## lm(formula = grp.mean ~ Female.condition, data = Condgroup_meanJW)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.6593 -0.2314 -0.0086 0.2142 0.6982
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.09533 0.66405 6.167 1.03e-05 ***
## Female.condition -0.02193 0.15616 -0.140 0.89
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3545 on 17 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.001158, Adjusted R-squared: -0.0576
## F-statistic: 0.01971 on 1 and 17 DF, p-value: 0.89No maternal traits considerred for potential covariate.
2. Final model and model performance
Model formula
lmer.JW = lmer(Juvenile.weight..mg.~ Dev.treatment * Repro.treatment + (1|Sump/Female.ID), data = Juv)## boundary (singular) fit: see help('isSingular')Performance checks
performance::check_model(lmer.JW, check="homogeneity") 
performance::check_model(lmer.JW, check="outliers")
performance::check_model(lmer.JW, check="qq", detrend = FALSE) 
performance::check_model(lmer.JW, check="normality") 
performance::check_model(lmer.JW, check="linearity") 
performance::check_model(lmer.JW, check="pp_check", re_formula=NA)
hist(residuals(lmer.JW), col="darkgray")
outlierTest(lmer.JW)## rstudent unadjusted p-value Bonferroni p
## 330 -3.872169 0.00012708 0.049053Fish #330 is light but not from a measurement issue. Kept in the analysis.
3. Summary statistics … cited in Table S4
summary(lmer.JW)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Juvenile.weight..mg. ~ Dev.treatment * Repro.treatment + (1 |
## Sump/Female.ID)
## Data: Juv
##
## REML criterion at convergence: 146.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.7766 -0.6094 0.0523 0.7163 2.8237
##
## Random effects:
## Groups Name Variance Std.Dev.
## Female.ID:Sump (Intercept) 1.154e-01 3.398e-01
## Sump (Intercept) 1.212e-10 1.101e-05
## Residual 7.090e-02 2.663e-01
## Number of obs: 386, groups: Female.ID:Sump, 20; Sump, 6
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 3.9619 0.1411 16.0739 28.088
## Dev.treatmentWarm 0.1756 0.2092 16.0552 0.840
## Repro.treatmentWarm -0.1688 0.1993 16.0268 -0.847
## Dev.treatmentWarm:Repro.treatmentWarm 0.1695 0.3213 16.0026 0.528
## Pr(>|t|)
## (Intercept) 4.32e-15 ***
## Dev.treatmentWarm 0.414
## Repro.treatmentWarm 0.410
## Dev.treatmentWarm:Repro.treatmentWarm 0.605
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Dv.trW Rpr.tW
## Dv.trtmntWr -0.674
## Rpr.trtmntW -0.708 0.477
## Dv.trtW:R.W 0.439 -0.651 -0.620
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')The model is singular as very little variance is assigned to Sump.
4. anova … cited in Table S7
anova(lmer.JW)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Dev.treatment 0.186212 0.186212 1 16.003 2.6264 0.1246
## Repro.treatment 0.019394 0.019394 1 16.003 0.2735 0.6081
## Dev.treatment:Repro.treatment 0.019737 0.019737 1 16.003 0.2784 0.60505. Emmeans
JWEMM = (emmeans(lmer.JW, ~ Dev.treatment*Repro.treatment) %>% as.data.frame)
JWEMM## Dev.treatment Repro.treatment emmean SE df lower.CL upper.CL
## Control Control 3.961949 0.1769743 3.82 3.461471 4.462427
## Warm Control 4.137534 0.1582119 9.24 3.781026 4.494041
## Control Warm 3.793171 0.1512178 6.08 3.424404 4.161939
## Warm Warm 4.138273 0.2252733 10.16 3.637413 4.639133
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.956. Plot … Figure 4c
The sample size for Warm-Control is 1 lower than SL and YA as 1 fish was removed as an outlier due to a measurement error (a separate fish from #330 identified above).
sample_size_Juv_NoOutlier <- c(94, 112, 60, 120)
fig.4c <- ggplot(JWEMM, aes(x = Repro.treatment, y = emmean, colour = Dev.treatment)) +
geom_jitter(data = Juv,
aes(x = Repro.treatment, y = Juvenile.weight..mg., 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, na.rm = TRUE) +
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 = 4, colour = "black", stroke = 0.5) +
geom_text(aes(label = sample_size_Juv_NoOutlier, 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 = "Weight (mg)") +
scale_fill_manual(values = c("Control" = "blue", "Warm" = "orange")) +
scale_colour_manual(values = c("Control" = "blue", "Warm" = "orange")) +
scale_y_continuous(breaks = c(3.0, 3.5, 4.0, 4.5, 5.0)) +
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.4c
Yolk area
1. Maternal effects and covariate exploration
SLgroup_meanJYa = ddply(Juv, "Female.SL..cm.", summarise, grp.mean=mean(Juvenile.yolk.area..mm.2.))
Wgroup_meanJYa = ddply(Juv, "Female.weight..g.", summarise, grp.mean=mean(Juvenile.yolk.area..mm.2.))
Condgroup_meanJYa = ddply(Juv, "Female.condition", summarise, grp.mean=mean(Juvenile.yolk.area..mm.2.))
Egggroup_meanYa = ddply(Juv, "Avg.egg.area", summarise, grp.mean=mean(Juvenile.yolk.area..mm.2.))Correlation with Female SL
lm.1 = lm(grp.mean ~ Female.SL..cm., data=SLgroup_meanJYa)
summary(lm.1)##
## Call:
## lm(formula = grp.mean ~ Female.SL..cm., data = SLgroup_meanJYa)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.33393 -0.14351 0.01688 0.10027 0.34179
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.57726 0.79116 3.258 0.00437 **
## Female.SL..cm. -0.13420 0.08892 -1.509 0.14859
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.18 on 18 degrees of freedom
## Multiple R-squared: 0.1123, Adjusted R-squared: 0.06301
## F-statistic: 2.278 on 1 and 18 DF, p-value: 0.1486Correlation with Female W
lm.2 = lm(grp.mean ~ Female.weight..g., data=Wgroup_meanJYa)
summary(lm.2)##
## Call:
## lm(formula = grp.mean ~ Female.weight..g., data = Wgroup_meanJYa)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.34399 -0.16740 0.03817 0.14471 0.33017
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.496589 0.282027 5.307 4.81e-05 ***
## Female.weight..g. -0.003779 0.009421 -0.401 0.693
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1902 on 18 degrees of freedom
## Multiple R-squared: 0.008859, Adjusted R-squared: -0.0462
## F-statistic: 0.1609 on 1 and 18 DF, p-value: 0.6931Correlation with Female FK
lm.3 = lm(grp.mean ~ Female.condition, data=Condgroup_meanJYa)
summary(lm.3)##
## Call:
## lm(formula = grp.mean ~ Female.condition, data = Condgroup_meanJYa)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.36281 -0.06290 -0.01345 0.13687 0.28644
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.83418 0.33167 2.515 0.0216 *
## Female.condition 0.13068 0.07815 1.672 0.1118
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1777 on 18 degrees of freedom
## Multiple R-squared: 0.1344, Adjusted R-squared: 0.08636
## F-statistic: 2.796 on 1 and 18 DF, p-value: 0.1118No maternal traits to be considered as potential covariate.
2. Final model and model performance
Model formula
lmer.JYA = lmer(Juvenile.yolk.area..mm.2.~ Dev.treatment * Repro.treatment + (1|Sump/Female.ID), data = Juv)## boundary (singular) fit: see help('isSingular')Performance checks
performance::check_model(lmer.JYA, check="homogeneity") 
performance::check_model(lmer.JYA, check="outliers")
performance::check_model(lmer.JYA, check="qq", detrend = FALSE) 
performance::check_model(lmer.JYA, check="normality") 
performance::check_model(lmer.JYA, check="linearity") 
performance::check_model(lmer.JYA, check="pp_check")
hist(residuals(lmer.JYA), col="darkgray")
outlierTest(lmer.JYA)## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
## rstudent unadjusted p-value Bonferroni p
## 166 3.762563 0.0001948 0.0753863. Summary statistics … cited in Table S7
summary(lmer.JYA)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Juvenile.yolk.area..mm.2. ~ Dev.treatment * Repro.treatment +
## (1 | Sump/Female.ID)
## Data: Juv
##
## REML criterion at convergence: -284.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.7663 -0.7306 -0.0595 0.6898 3.6693
##
## Random effects:
## Groups Name Variance Std.Dev.
## Female.ID:Sump (Intercept) 3.599e-02 1.897e-01
## Sump (Intercept) 1.247e-11 3.531e-06
## Residual 2.299e-02 1.516e-01
## Number of obs: 387, groups: Female.ID:Sump, 20; Sump, 6
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 1.410625 0.078806 16.081892 17.900
## Dev.treatmentWarm 0.044268 0.116852 16.062548 0.379
## Repro.treatmentWarm -0.099300 0.111355 16.028668 -0.892
## Dev.treatmentWarm:Repro.treatmentWarm 0.006192 0.179491 16.006260 0.034
## Pr(>|t|)
## (Intercept) 4.81e-12 ***
## Dev.treatmentWarm 0.710
## Repro.treatmentWarm 0.386
## Dev.treatmentWarm:Repro.treatmentWarm 0.973
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Dv.trW Rpr.tW
## Dv.trtmntWr -0.674
## Rpr.trtmntW -0.708 0.477
## Dv.trtW:R.W 0.439 -0.651 -0.620
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')4. anova … cited in Table S4
anova(lmer.JYA)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Dev.treatment 0.0064035 0.0064035 1 16.006 0.2785 0.6049
## Repro.treatment 0.0264186 0.0264186 1 16.006 1.1491 0.2996
## Dev.treatment:Repro.treatment 0.0000274 0.0000274 1 16.006 0.0012 0.97295. Emmeans
JYEMM = (emmeans(lmer.JYA, ~ Dev.treatment*Repro.treatment) %>% as.data.frame)
JYEMM## Dev.treatment Repro.treatment emmean SE df lower.CL upper.CL
## Control Control 1.410625 0.09886771 3.82 1.131088 1.690162
## Warm Control 1.454893 0.08839064 9.24 1.255726 1.654060
## Control Warm 1.311325 0.08446638 6.08 1.105326 1.517324
## Warm Warm 1.361785 0.12585658 10.16 1.081964 1.641605
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95emmeans(lmer.JYA, pairwise ~ Dev.treatment * Repro.treatment, adjust="tukey")## $emmeans
## Dev.treatment Repro.treatment emmean SE df lower.CL upper.CL
## Control Control 1.41 0.0989 3.82 1.13 1.69
## Warm Control 1.45 0.0884 9.24 1.26 1.65
## Control Warm 1.31 0.0845 6.08 1.11 1.52
## Warm Warm 1.36 0.1259 10.16 1.08 1.64
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## Control Control - Warm Control -0.0443 0.128 14.87 -0.346 0.9852
## Control Control - Control Warm 0.0993 0.130 4.79 0.764 0.8675
## Control Control - Warm Warm 0.0488 0.160 7.24 0.305 0.9893
## Warm Control - Control Warm 0.1436 0.122 7.63 1.174 0.6587
## Warm Control - Warm Warm 0.0931 0.154 9.81 0.605 0.9280
## Control Warm - Warm Warm -0.0505 0.143 14.73 -0.353 0.9843
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 4 estimates6. Plot … Figure 4d
fig.4d <- ggplot(JYEMM, aes(x = Repro.treatment, y = emmean, colour = Dev.treatment)) +
geom_jitter(data = Juv,
aes(x = Repro.treatment, y = Juvenile.yolk.area..mm.2., 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 = 4, colour = "black", stroke = 0.5) +
geom_text(aes(label = sample_size_Juv, 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 = "Yolk area (mm²)") +
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.4d
Condition: W for a givel SL
SL and W both did not have maternal traits assigned as a covariate. The model for Condition is run without any covariates.
1. Final model and model performance
Model formula
lmer.JCondlog = lmer(log.W ~ Dev.treatment * Repro.treatment + log.SL + (1|Sump/Female.ID) , data = Juv) Performance checks
performance::check_model(lmer.JCondlog, check="homogeneity") 
performance::check_model(lmer.JCondlog, check="outliers")
performance::check_model(lmer.JCondlog, check="qq", detrend = FALSE) 
performance::check_model(lmer.JCondlog, check="normality") 
performance::check_model(lmer.JCondlog, check="linearity") 
performance::check_model(lmer.JCondlog, check="pp_check", re_formula=NA)
hist(residuals(lmer.JCondlog), col="darkgray")
outlierTest(lmer.JCondlog)## rstudent unadjusted p-value Bonferroni p
## 69 -5.797071 1.4277e-08 5.5109e-06
## 65 5.111673 5.0879e-07 1.9639e-04
## 238 4.687180 3.8761e-06 1.4962e-03Similar to W, the fish identified as outliers (fish 65, 69 and 238) were produced by pairs that generally produced heavier hatchlings for a given length and were retained in the analysis.
2. Summary statistics … cited in Table S7
summary(lmer.JCondlog)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## log.W ~ Dev.treatment * Repro.treatment + log.SL + (1 | Sump/Female.ID)
## Data: Juv
##
## REML criterion at convergence: -1978.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.6528 -0.4541 -0.0351 0.4478 4.9686
##
## Random effects:
## Groups Name Variance Std.Dev.
## Female.ID:Sump (Intercept) 2.304e-03 4.800e-02
## Sump (Intercept) 2.636e-11 5.134e-06
## Residual 2.510e-04 1.584e-02
## Number of obs: 386, groups: Female.ID:Sump, 20; Sump, 6
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) -0.33567 0.03687 151.41935 -9.104
## Dev.treatmentWarm 0.01408 0.02915 15.93407 0.483
## Repro.treatmentWarm -0.01613 0.02779 15.92855 -0.580
## log.SL 1.29726 0.04344 369.02782 29.866
## Dev.treatmentWarm:Repro.treatmentWarm 0.03545 0.04482 15.92826 0.791
## Pr(>|t|)
## (Intercept) 4.76e-16 ***
## Dev.treatmentWarm 0.636
## Repro.treatmentWarm 0.570
## log.SL < 2e-16 ***
## Dev.treatmentWarm:Repro.treatmentWarm 0.441
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Dv.trW Rpr.tW log.SL
## Dv.trtmntWr -0.355
## Rpr.trtmntW -0.380 0.477
## log.SL -0.846 -0.005 0.004
## Dv.trtW:R.W 0.224 -0.651 -0.620 0.0123. anova … cited in Table S4
anova(lmer.JCondlog)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Dev.treatment 0.000506 0.000506 1 15.92 2.0145 0.1751
## Repro.treatment 0.000001 0.000001 1 15.93 0.0051 0.9440
## log.SL 0.223897 0.223897 1 369.03 891.9794 <2e-16
## Dev.treatment:Repro.treatment 0.000157 0.000157 1 15.93 0.6257 0.4406
##
## Dev.treatment
## Repro.treatment
## log.SL ***
## Dev.treatment:Repro.treatment
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 14. emmeans
Since SL and W used to estimate Condition was log10() transformed, Emmeans and the 95% confidence intervals are back-transformed with 10^(). SE is removed from the table as they are not in the same scale as the response.
JCondEMM <- emmeans(lmer.JCondlog, ~ Dev.treatment * Repro.treatment) %>%
as.data.frame() %>%
mutate(emmean = 10^(emmean),
lower.CL = 10^(lower.CL),
upper.CL = 10^(upper.CL)) %>%
select(-SE)
JCondEMM## Dev.treatment Repro.treatment emmean df lower.CL upper.CL
## 1 Control Control 3.922264 3.791470 3.337845 4.609008
## 2 Warm Control 4.051495 9.208026 3.613109 4.543073
## 3 Control Warm 3.779287 6.088331 3.357151 4.254504
## 4 Warm Warm 4.235826 10.179629 3.606259 4.9753005. Plot … Figure 4b
Jiterred points show fitted values of Condition estimated for each
individual fish.
The sample size in Control-Warm is 1 lower than SL and yolk area due to
a measurement error in W for one fish.
Juv$Cond_fitted <- 10^(predict(lmer.JCondlog, newdata = Juv, re.form = NULL))fig.4b <- ggplot(JCondEMM, aes(x = Repro.treatment, y = emmean, colour = Dev.treatment)) +
geom_jitter(data = Juv,
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, na.rm = TRUE) +
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 = 4, colour = "black", stroke = 0.5) +
geom_text(aes(label = sample_size_Juv_NoOutlier, 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 = "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.4b
How egg size infleunces thermal effects on newly hatched traits
Egg size is expected to influence juvenile quality at hatching.
Below, the correlation between egg size and newly hatched traits are
first investigated using linear regression. The impact of egg size on
thermal effects are then evaluated for traits correlated with egg
size.
Eggareagroup_meanJSL = ddply(Juv, "Avg.egg.area", summarise, grp.mean=mean(Juvenile.SL))
Egggroup_meanJW = ddply(Juv, "Avg.egg.area", summarise, grp.mean=mean(Juvenile.weight..mg.))
Egggroup_meanYa = ddply(Juv, "Avg.egg.area", summarise, grp.mean=mean(Juvenile.yolk.area..mm.2.))Correlation of newly hatched traits with egg size
SL
lm.SL = lm(grp.mean ~ Avg.egg.area, data=Eggareagroup_meanJSL)
summary(lm.SL)##
## Call:
## lm(formula = grp.mean ~ Avg.egg.area, data = Eggareagroup_meanJSL)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.3174 -0.1799 0.0255 0.1260 0.3864
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.81657 0.45400 8.407 1.2e-07 ***
## Avg.egg.area 0.29707 0.09597 3.095 0.00624 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2024 on 18 degrees of freedom
## Multiple R-squared: 0.3474, Adjusted R-squared: 0.3111
## F-statistic: 9.581 on 1 and 18 DF, p-value: 0.006243W
lm.W = lm(grp.mean ~ Avg.egg.area, data=Egggroup_meanJW)
summary(lm.W)##
## Call:
## lm(formula = grp.mean ~ Avg.egg.area, data = Egggroup_meanJW)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.69584 -0.22947 -0.01227 0.22305 0.70013
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.30416 0.82207 5.236 6.71e-05 ***
## Avg.egg.area -0.06367 0.17283 -0.368 0.717
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3533 on 17 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.00792, Adjusted R-squared: -0.05044
## F-statistic: 0.1357 on 1 and 17 DF, p-value: 0.7171Yolk area
lm.YA = lm(grp.mean ~ Avg.egg.area, data=Egggroup_meanYa)
summary(lm.YA)##
## Call:
## lm(formula = grp.mean ~ Avg.egg.area, data = Egggroup_meanYa)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.264078 -0.117863 0.009804 0.068614 0.254205
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.30291 0.34340 0.882 0.38935
## Avg.egg.area 0.22985 0.07259 3.166 0.00535 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1531 on 18 degrees of freedom
## Multiple R-squared: 0.3577, Adjusted R-squared: 0.322
## F-statistic: 10.02 on 1 and 18 DF, p-value: 0.005345Juvenile SL and YA correlated with egg size.
How egg size influences thermal impacts on these traits is explored
below.
Does egg size influence thermal effects on juvenile SL?
lmer.JSLb = lmer(Juvenile.SL~ Dev.treatment * Repro.treatment + Avg.egg.area + (1 |Sump/Female.ID), data = Juv)## boundary (singular) fit: see help('isSingular')Summary statistics
lmer.JSLb = lmer(Juvenile.SL~ Dev.treatment * Repro.treatment + Avg.egg.area + (1 |Sump/Female.ID), data = Juv)## boundary (singular) fit: see help('isSingular')summary(lmer.JSLb)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Juvenile.SL ~ Dev.treatment * Repro.treatment + Avg.egg.area +
## (1 | Sump/Female.ID)
## Data: Juv
##
## REML criterion at convergence: 7.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3257 -0.5655 0.0151 0.6994 2.7257
##
## Random effects:
## Groups Name Variance Std.Dev.
## Female.ID:Sump (Intercept) 0.04104 0.2026
## Sump (Intercept) 0.00000 0.0000
## Residual 0.05042 0.2245
## Number of obs: 387, groups: Female.ID:Sump, 20; Sump, 6
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 3.43396 0.59502 14.99429 5.771
## Dev.treatmentWarm 0.06849 0.12702 14.98707 0.539
## Repro.treatmentWarm 0.17494 0.13890 15.01384 1.259
## Avg.egg.area 0.36232 0.11839 14.95298 3.060
## Dev.treatmentWarm:Repro.treatmentWarm -0.20825 0.19556 14.89024 -1.065
## Pr(>|t|)
## (Intercept) 3.7e-05 ***
## Dev.treatmentWarm 0.59764
## Repro.treatmentWarm 0.22708
## Avg.egg.area 0.00796 **
## Dev.treatmentWarm:Repro.treatmentWarm 0.30388
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Dv.trW Rpr.tW Avg.g.
## Dv.trtmntWr -0.167
## Rpr.trtmntW -0.578 0.449
## Avg.egg.are -0.990 0.071 0.495
## Dv.trtW:R.W 0.167 -0.654 -0.588 -0.105
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')Singular as no varriance is assigned to Sump.
Comparing anova
Final model (without average egg size)
anova(lmer.JSL)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Dev.treatment 0.006708 0.006708 1 14.3039 0.1331 0.7206
## Repro.treatment 0.046344 0.046344 1 2.5837 0.9194 0.4185
## Dev.treatment:Repro.treatment 0.019615 0.019615 1 14.3039 0.3891 0.5426Model with average egg size
lmer.JSLb = lmer(Juvenile.SL~ Dev.treatment * Repro.treatment + Avg.egg.area + (1 |Sump/Female.ID), data = Juv)## boundary (singular) fit: see help('isSingular')anova(lmer.JSLb)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Dev.treatment 0.00677 0.00677 1 14.881 0.1343 0.719213
## Repro.treatment 0.01964 0.01964 1 14.979 0.3896 0.541905
## Avg.egg.area 0.47224 0.47224 1 14.953 9.3664 0.007956 **
## Dev.treatment:Repro.treatment 0.05718 0.05718 1 14.890 1.1340 0.303878
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Avg egg size significantly influences juvenile SL (larger eggs lead to longer hatchlings), but the overall patterns of thermal effects remains the same.
Does egg size influence thermal effects on yolk area?
lmer.JYAb = lmer(Juvenile.yolk.area..mm.2.~ Dev.treatment * Repro.treatment + Avg.egg.area + (1|Sump/Female.ID), data = Juv)Summary statistics
lmer.JYAb = lmer(Juvenile.yolk.area..mm.2.~ Dev.treatment * Repro.treatment + Avg.egg.area + (1|Sump/Female.ID), data = Juv)
summary(lmer.JYAb)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Juvenile.yolk.area..mm.2. ~ Dev.treatment * Repro.treatment +
## Avg.egg.area + (1 | Sump/Female.ID)
## Data: Juv
##
## REML criterion at convergence: -288.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.7656 -0.7309 -0.0613 0.6728 3.6766
##
## Random effects:
## Groups Name Variance Std.Dev.
## Female.ID:Sump (Intercept) 0.02452 0.15658
## Sump (Intercept) 0.00191 0.04371
## Residual 0.02300 0.15164
## Number of obs: 387, groups: Female.ID:Sump, 20; Sump, 6
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 0.17707 0.46948 14.81364 0.377
## Dev.treatmentWarm 0.06769 0.09850 12.06794 0.687
## Repro.treatmentWarm 0.05471 0.11513 3.84474 0.475
## Avg.egg.area 0.24647 0.09284 14.74834 2.655
## Dev.treatmentWarm:Repro.treatmentWarm -0.05112 0.15160 11.93265 -0.337
## Pr(>|t|)
## (Intercept) 0.7114
## Dev.treatmentWarm 0.5049
## Repro.treatmentWarm 0.6603
## Avg.egg.area 0.0182 *
## Dev.treatmentWarm:Repro.treatmentWarm 0.7418
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Dv.trW Rpr.tW Avg.g.
## Dv.trtmntWr -0.180
## Rpr.trtmntW -0.572 0.437
## Avg.egg.are -0.988 0.084 0.481
## Dv.trtW:R.W 0.165 -0.654 -0.552 -0.103Singular as very little variance assigned to Sump.
Comparing anova
Final model (without egg size included)
anova(lmer.JYA)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Dev.treatment 0.0064035 0.0064035 1 16.006 0.2785 0.6049
## Repro.treatment 0.0264186 0.0264186 1 16.006 1.1491 0.2996
## Dev.treatment:Repro.treatment 0.0000274 0.0000274 1 16.006 0.0012 0.9729Model with egg size included
lmer.JYAb = lmer(Juvenile.yolk.area..mm.2.~ Dev.treatment * Repro.treatment + Avg.egg.area + (1|Sump/Female.ID), data = Juv)
anova(lmer.JYAb)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Dev.treatment 0.007180 0.007180 1 12.080 0.3122 0.58653
## Repro.treatment 0.002088 0.002088 1 1.956 0.0908 0.79216
## Avg.egg.area 0.162065 0.162065 1 14.748 7.0477 0.01822 *
## Dev.treatment:Repro.treatment 0.002615 0.002615 1 11.933 0.1137 0.74181
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Average egg size significantly influences newly hatched yolk area (larger eggs lead to newly hatched fish having more yolk reserve), but does not change the overall pattern of thermal impacts.
Combining plots for Figure 4
fig.4_legend <- fig.4a + theme(legend.position = "top")
legend <- get_legend(fig.4_legend)
fig.Juv <- ggarrange(
fig.4a, fig.4b, fig.4c, fig.4d,
labels = c("a", "b", "c", "d"),
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.Juv_with_legend <- ggarrange(
legend, fig.Juv,
ncol = 1, heights = c(0.1, 0.9)
)
fig.Juv_with_legend <- annotate_figure(
fig.Juv_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 4.pdf",
plot = fig.Juv_with_legend,
width = 169,
height = 180,
units = "mm",
dpi = 600
)