Thermal history influences reproduction and offspring quality in a coral reef fish: Reproductive performance
Y.Yasutake
March 2022, updated January 2025
Rmarkdown
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 loading
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$Family.combo = factor(Parental$Family.combo)
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 : Factor w/ 1 level "AC": 1 1 1 1 1 1 1 1 1 1 ...
## $ 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 ...Clutches per pair
Clutchpp <- read_excel("ThermalHistoryCoralReefFish.xlsx",
sheet = "Clutches per pair",
trim_ws = TRUE) %>%
setNames(make.names(names(.), unique = TRUE))
Clutchpp$Sump = factor(Clutchpp$Sump)
Clutchpp$Tank.ID = factor(Clutchpp$Tank.ID)
Clutchpp$Treatment = factor(Clutchpp$Treatment)
Clutchpp$Dev.treatment = factor(Clutchpp$Dev.treatment)
Clutchpp$Repro.treatment = factor(Clutchpp$Repro.treatment)
Clutchpp$Female.family = factor(Clutchpp$Male.family)
Clutchpp$Male.family = factor(Clutchpp$Male.family)
Clutchpp$Male.family = factor(Clutchpp$Male.family)
Clutchpp$Family.combo = factor(Clutchpp$Family.combo)
str(Clutchpp)## tibble [20 × 18] (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 : chr [1:20] "Control-Control-1" "Control-Control-2" "Control-Control-3" "Control-Control-4" ...
## $ Family.combo : Factor w/ 11 levels "AC","AD","AE",..: 4 2 10 11 1 9 3 2 8 7 ...
## $ Female.ID : chr [1:20] "ACC4" "ACC3" "ECCb" "ECC" ...
## $ Female.family : Factor w/ 5 levels "A","C","D","E",..: 5 3 1 5 2 1 4 3 4 3 ...
## $ Female.SL..cm. : num [1:20] 9.14 8.91 9.27 7.8 8.78 ...
## $ Female.W..g. : num [1:20] 31.3 35.9 33 26.6 32.6 ...
## $ Female.condition : 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 ...
## $ Male.SL..cm. : num [1:20] 8.65 9.65 9.45 7.74 8.6 ...
## $ Male.weight..g. : num [1:20] 22.6 41.1 36.9 19.7 23.6 ...
## $ Male.condition : num [1:20] 3.5 4.58 4.37 4.25 3.71 ...
## $ Number.of.clutches.per.pair: num [1:20] 3 3 2 3 3 2 3 3 3 3 ...Clutch size
ClutchSz <- read_excel("ThermalHistoryCoralReefFish.xlsx",
sheet = "Clutch size",
trim_ws = TRUE) %>%
setNames(make.names(names(.), unique = TRUE))
ClutchSz$Sump = factor(ClutchSz$Sump)
ClutchSz$Tank.ID = factor(ClutchSz$Tank.ID)
ClutchSz$Treatment = factor(ClutchSz$Treatment)
ClutchSz$Dev.treatment = factor(ClutchSz$Dev.treatment)
ClutchSz$Repro.treatment = factor(ClutchSz$Repro.treatment)
ClutchSz$Male.family = factor(ClutchSz$Male.family)
ClutchSz$Female.family = factor(ClutchSz$Female.family)
ClutchSz$Family.combo = factor(ClutchSz$Family.combo)
str(ClutchSz)## tibble [20 × 18] (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 : chr [1:20] "Control-Control-1" "Control-Control-2" "Control-Control-3" "Control-Control-4" ...
## $ Family.combo : Factor w/ 11 levels "AC","AD","AE",..: 4 2 10 11 1 9 3 2 8 7 ...
## $ Female.ID : chr [1:20] "ACC4" "ACC3" "ECCb" "ECC" ...
## $ Female.family : Factor w/ 5 levels "A","B","C","D",..: 1 1 5 5 1 4 1 1 3 3 ...
## $ Female.SL..cm. : num [1:20] 9.14 8.91 9.27 7.8 8.78 ...
## $ Female.Weight..g.: num [1:20] 31.3 35.9 33 26.6 32.6 ...
## $ Female.condition : 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 ...
## $ Male.SL..cm. : num [1:20] 8.65 9.65 9.45 7.74 8.6 ...
## $ Male.Weight..g. : num [1:20] 22.6 41.1 36.9 19.7 23.6 ...
## $ Male.condition : num [1:20] 3.5 4.58 4.37 4.25 3.71 ...
## $ Clutch.size..1st.: num [1:20] 198 324 309 123 349 302 216 323 364 278 ...Total number of eggs per pair
Totals <- read_excel("ThermalHistoryCoralReefFish.xlsx",
sheet = "Total eggs per pair",
trim_ws = TRUE) %>%
setNames(make.names(names(.), unique = TRUE))
Totals$Sump = factor(Totals$Sump)
Totals$Tank.ID = factor(Totals$Tank.ID)
Totals$Treatment = factor(Totals$Treatment)
Totals$Dev.treatment = factor(Totals$Dev.treatment)
Totals$Repro.treatment = factor(Totals$Repro.treatment)
Totals$Male.family = factor(Totals$Male.family)
Totals$Female.family = factor(Totals$Female.family)
Totals$Family.combo = factor(Totals$Family.combo)
str(Totals)## tibble [20 × 18] (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 : chr [1:20] "Control-Control-1" "Control-Control-2" "Control-Control-3" "Control-Control-4" ...
## $ Family.combo : Factor w/ 11 levels "AC","AD","AE",..: 4 2 10 11 1 9 3 2 8 7 ...
## $ Female.ID : chr [1:20] "ACC4" "ACC3" "ECCb" "ECC" ...
## $ Female.family : Factor w/ 5 levels "A","B","C","D",..: 1 1 5 5 1 4 1 1 3 3 ...
## $ Female.SL..cm. : num [1:20] 9.14 8.91 9.27 7.8 8.78 ...
## $ Female.W..g. : num [1:20] 31.3 35.9 33 26.6 32.6 ...
## $ Female.condition : 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 ...
## $ Male.SL..cm. : num [1:20] 8.65 9.65 9.45 7.74 8.6 ...
## $ Male.weight..g. : num [1:20] 22.6 41.1 36.9 19.7 23.6 ...
## $ Male.condition : 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 ...Egg size
EggArea <- read_excel("ThermalHistoryCoralReefFish.xlsx",
sheet = "Egg size",
trim_ws = TRUE) %>%
setNames(make.names(names(.), unique = TRUE))
EggArea$Sump = factor(EggArea$Sump)
EggArea$Tank.ID = factor(EggArea$Tank.ID)
EggArea$Treatment = factor(EggArea$Treatment)
EggArea$Dev.treatment = factor(EggArea$Dev.treatment)
EggArea$Repro.treatment = factor(EggArea$Repro.treatment)
EggArea$Male.family = factor(EggArea$Male.family)
EggArea$Female.family = factor(EggArea$Female.family)
EggArea$Family.combo = factor(EggArea$Family.combo)
EggArea$Egg.replicate = factor(EggArea$Egg.replicate)
EggArea$Clutch.replicate = factor(EggArea$Clutch.replicate)
str(EggArea)## tibble [199 × 19] (S3: tbl_df/tbl/data.frame)
## $ Tank.ID : Factor w/ 20 levels "9","25","42",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ 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 ...
## $ Egg.replicate : Factor w/ 199 levels "Control_Control_1_1",..: 170 172 173 174 175 176 177 178 179 171 ...
## $ Female.ID : chr [1:199] "AHH" "AHH" "AHH" "AHH" ...
## $ Male.ID : chr [1:199] "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 : Factor w/ 11 levels "AC","AD","AE",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Male.SL..cm. : num [1:199] 9.48 9.48 9.48 9.48 9.48 ...
## $ Male.weight..g. : num [1:199] 31.9 31.9 31.9 31.9 31.9 31.9 31.9 31.9 31.9 31.9 ...
## $ Male.condition : num [1:199] 3.75 3.75 3.75 3.75 3.75 ...
## $ Female.SL..cm. : num [1:199] 9.07 9.07 9.07 9.07 9.07 ...
## $ Female.weight..g.: num [1:199] 33.2 33.2 33.2 33.2 33.2 33.2 33.2 33.2 33.2 33.2 ...
## $ Female.condition : num [1:199] 4.44 4.44 4.44 4.44 4.44 ...
## $ Egg.area..mm2. : num [1:199] 4.99 4.81 5.17 5.28 4.96 ...Individual trait models
For each reproductive performance trait, the information is presented
as follows:
1. Error structure assessment: A dispersion is used to
test whether a Poisson or Negative Binomial distribution is appropriate
for count data (clutches per pair, clutch size, and total number of eggs
per pair). For non-count data, a Gaussian distribution is used.
2. Maternal effects and covariate exploration:
A linear regression is performed between each reproductive trait and
maternal physical variables.
A maternal variable is included in the final model if it improves model
fit based on AIC values.
3. Final model and model performance:
4. Summary statistics:
Summary statistics are cited in Table S6 for each trait.
5. anova results:
Results from anova are cited in Table S3 for each trait.
6. Estimated marginal means (emmeans):
Emmeans and 95% confidence intervals are provided for each trait. These
values are used to generate the corresponding plots.
7. Trait plots:
Plots for each traits are provided, which are combined at the end of the
markdown to create Figure 4.
Clutches per pair
1. Error structure assessment
poisson.clutchpp <- glmmTMB(Number.of.clutches.per.pair ~ Dev.treatment * Repro.treatment, data=Clutchpp, family="poisson")negb.clutchpp <- glmmTMB(Number.of.clutches.per.pair ~ Dev.treatment * Repro.treatment, data=Clutchpp, family="nbinom1")## Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
## problem; non-positive-definite Hessian matrix. See vignette('troubleshooting')Negative binomial has issues with convergence. The Poisson model is chosen.
res.poisson.clutchpp <- simulateResiduals(poisson.clutchpp, plot=T)
Poisson has minor issues with dispersion, but this is preferable over
the convergence issue with the negative binomial.
Proceeding with Poisson keeping the dispersion issue in mind.
2. Exploring maternal effects
Correlation with Female SL
lm.1 = lm(Number.of.clutches.per.pair ~ SL...Female..cm., data=Parental)
summary(lm.1)##
## Call:
## lm(formula = Number.of.clutches.per.pair ~ SL...Female..cm.,
## data = Parental)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.3807 -0.3522 0.1583 0.7089 0.7602
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.4829 3.6080 0.965 0.347
## SL...Female..cm. -0.1331 0.4055 -0.328 0.746
##
## Residual standard error: 0.8208 on 18 degrees of freedom
## Multiple R-squared: 0.005951, Adjusted R-squared: -0.04927
## F-statistic: 0.1078 on 1 and 18 DF, p-value: 0.7465Correlation with Female W
lm.2 = lm(Number.of.clutches.per.pair ~ Weight...Female..g., data=Parental)
summary(lm.2)##
## Call:
## lm(formula = Number.of.clutches.per.pair ~ Weight...Female..g.,
## data = Parental)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.5414 -0.3987 0.1905 0.6905 0.8218
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.55768 1.20792 1.290 0.214
## Weight...Female..g. 0.02508 0.04035 0.622 0.542
##
## Residual standard error: 0.8146 on 18 degrees of freedom
## Multiple R-squared: 0.02102, Adjusted R-squared: -0.03337
## F-statistic: 0.3865 on 1 and 18 DF, p-value: 0.542Correlation with Female FK
lm.3 = lm(Number.of.clutches.per.pair ~ Condition...Female, data=Parental)
summary(lm.3)##
## Call:
## lm(formula = Number.of.clutches.per.pair ~ Condition...Female,
## data = Parental)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.22395 -0.46314 -0.02657 0.63694 1.01189
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.08599 1.42808 -0.060 0.953
## Condition...Female 0.56630 0.33650 1.683 0.110
##
## Residual standard error: 0.7653 on 18 degrees of freedom
## Multiple R-squared: 0.136, Adjusted R-squared: 0.08795
## F-statistic: 2.832 on 1 and 18 DF, p-value: 0.1097No maternal traits to be considered for a covariate.
3. Final Model
Model formula
poisson.clutchpp <- glmmTMB(Number.of.clutches.per.pair ~ Dev.treatment * Repro.treatment, data=Clutchpp, family="poisson")Performance checks
performance::check_zeroinflation(poisson.clutchpp)## Model has no observed zeros in the response variable.## NULLperformance::check_model(poisson.clutchpp, residual_type = "response")## `check_outliers()` does not yet support models of class `glmmTMB`.
hist(residuals(poisson.clutchpp), col="darkgray")
4. Summary statistics … cited in Table S3
summary(poisson.clutchpp)## Family: poisson ( log )
## Formula: Number.of.clutches.per.pair ~ Dev.treatment * Repro.treatment
## Data: Clutchpp
##
## AIC BIC logLik deviance df.resid
## 66.3 70.3 -29.2 58.3 16
##
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.98083 0.25000 3.923 8.73e-05 ***
## Dev.treatmentWarm -0.28768 0.40311 -0.714 0.475
## Repro.treatmentWarm -0.06454 0.35940 -0.180 0.857
## Dev.treatmentWarm:Repro.treatmentWarm -0.11778 0.65511 -0.180 0.857
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 15. Anova … cited in Table S6
Anova(poisson.clutchpp)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Number.of.clutches.per.pair
## Chisq Df Pr(>Chisq)
## Dev.treatment 1.0935 1 0.2957
## Repro.treatment 0.1107 1 0.7393
## Dev.treatment:Repro.treatment 0.0323 1 0.85736. Emmeans
ClutchppEMM = (emmeans(poisson.clutchpp, ~ Dev.treatment * Repro.treatment, type = "response") %>% as.data.frame)
ClutchppEMM## Dev.treatment Repro.treatment rate SE df asymp.LCL asymp.UCL
## Control Control 2.666667 0.6666667 Inf 1.6336858 4.352805
## Warm Control 2.000000 0.6324554 Inf 1.0761092 3.717093
## Control Warm 2.500000 0.6454972 Inf 1.5071648 4.146859
## Warm Warm 1.666669 0.7453565 Inf 0.6937145 4.004220
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scalePlot …. Figure 3a
For all trait plots, post-maturation treatments are given as the
x-axis, and 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 from each breeding pair shown
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(6, 5, 6, 3)
fig.3a <- ggplot(ClutchppEMM, aes(x = Repro.treatment, y = rate, colour = Dev.treatment)) +
geom_jitter(data = Clutchpp,
aes(x = Repro.treatment, y = Number.of.clutches.per.pair, colour = Dev.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 = asymp.LCL, ymax = asymp.UCL, 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_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 = "Cluthces per pair") +
scale_y_continuous(breaks = c(1.0, 2.0, 3.0)) +
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.3a
Clutch size
1. Error structure assessment
poisson.Csz <- glmmTMB(Clutch.size..1st. ~ Dev.treatment * Repro.treatment, data=ClutchSz, family="poisson")
negb.Csz <- glmmTMB(Clutch.size..1st. ~ Dev.treatment * Repro.treatment, data=ClutchSz, family="nbinom1")res.poisson.Csz<-simulateResiduals(poisson.Csz, plot=T)
res.negb.Csz<-simulateResiduals(negb.Csz, plot=T)
Poisson has issues with dispersion. Choosing the negative binomial model.
2. Maternal effects
Correlation to Female SL
lm.1 = lm(Clutch.size..1st. ~ SL...Female..cm., data=Parental)
summary(lm.1)##
## Call:
## lm(formula = Clutch.size..1st. ~ SL...Female..cm., data = Parental)
##
## Residuals:
## Min 1Q Median 3Q Max
## -104.216 -36.695 4.832 42.940 99.791
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -456.78 249.58 -1.830 0.0838 .
## SL...Female..cm. 80.45 28.05 2.868 0.0102 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 56.78 on 18 degrees of freedom
## Multiple R-squared: 0.3136, Adjusted R-squared: 0.2755
## F-statistic: 8.224 on 1 and 18 DF, p-value: 0.01023Correlation to Female W
lm.2 = lm(Clutch.size..1st. ~ Weight...Female..g., data=Parental)
summary(lm.2)##
## Call:
## lm(formula = Clutch.size..1st. ~ Weight...Female..g., data = Parental)
##
## Residuals:
## Min 1Q Median 3Q Max
## -114.44 -27.55 11.96 34.57 102.45
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 52.976 89.093 0.595 0.5595
## Weight...Female..g. 6.929 2.976 2.328 0.0317 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 60.08 on 18 degrees of freedom
## Multiple R-squared: 0.2315, Adjusted R-squared: 0.1888
## F-statistic: 5.422 on 1 and 18 DF, p-value: 0.03174Correlation to Female FK
lm.3 = lm(Clutch.size..1st. ~ Condition...Female, data=Parental)
summary(lm.3)##
## Call:
## lm(formula = Clutch.size..1st. ~ Condition...Female, data = Parental)
##
## Residuals:
## Min 1Q Median 3Q Max
## -114.47 -47.60 -10.67 61.95 99.98
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 320.16 127.04 2.520 0.0214 *
## Condition...Female -14.74 29.93 -0.492 0.6284
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 68.08 on 18 degrees of freedom
## Multiple R-squared: 0.01329, Adjusted R-squared: -0.04152
## F-statistic: 0.2425 on 1 and 18 DF, p-value: 0.6284Female SL has the strongest correlation. Using it as the potential covariate.
3. Covariate testing
negb.Csz <- glmmTMB(Clutch.size..1st. ~ Dev.treatment * Repro.treatment, data=ClutchSz, family="nbinom1")
negb.Csza <- glmmTMB(Clutch.size..1st. ~ Dev.treatment * Repro.treatment + Female.SL..cm., data=ClutchSz, family="nbinom1")AIC(negb.Csz, negb.Csza)## df AIC
## negb.Csz 5 231.3107
## negb.Csza 6 227.8299Maternal SL improves model fit, including in the final model.
4. Final model
Model formula
negb.Csza <- glmmTMB(Clutch.size..1st. ~ Dev.treatment * Repro.treatment + Female.SL..cm., data=ClutchSz, family="nbinom1")Performance check
performance::check_overdispersion(negb.Csza)## # Overdispersion test
##
## dispersion ratio = 1.055
## p-value = 0.744## No overdispersion detected.performance::check_model(negb.Csza, residual_type = "response")## `check_outliers()` does not yet support models of class `glmmTMB`.
hist(residuals(negb.Csza), col="darkgray")
5. Summary statistics … cited in Table S6
summary(negb.Csza)## Family: nbinom1 ( log )
## Formula:
## Clutch.size..1st. ~ Dev.treatment * Repro.treatment + Female.SL..cm.
## Data: ClutchSz
##
## AIC BIC logLik deviance df.resid
## 227.8 233.8 -107.9 215.8 14
##
##
## Dispersion parameter for nbinom1 family (): 10.5
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.143449 0.995119 3.159 0.00158 **
## Dev.treatmentWarm -0.123816 0.132724 -0.933 0.35088
## Repro.treatmentWarm 0.007446 0.123437 0.060 0.95190
## Female.SL..cm. 0.273266 0.111368 2.454 0.01414 *
## Dev.treatmentWarm:Repro.treatmentWarm 0.111384 0.204672 0.544 0.58630
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 16. Anova … cited in Table S3
Anova(negb.Csza)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Clutch.size..1st.
## Chisq Df Pr(>Chisq)
## Dev.treatment 0.5744 1 0.44852
## Repro.treatment 0.2607 1 0.60965
## Female.SL..cm. 6.0207 1 0.01414 *
## Dev.treatment:Repro.treatment 0.2962 1 0.58630
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 17. Emmeans
CszEMM <- emmeans(negb.Csza, ~ Dev.treatment * Repro.treatment, type = "response") %>%
as.data.frame()
CszEMM ## Dev.treatment Repro.treatment response SE df asymp.LCL asymp.UCL
## Control Control 262.8717 22.20804 Inf 222.7575 310.2097
## Warm Control 232.2584 23.80404 Inf 189.9906 283.9295
## Control Warm 264.8364 24.16885 Inf 221.4610 316.7072
## Warm Warm 261.5644 33.36982 Inf 203.6968 335.8715
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale8. Plot … Figure 3b
sample_sizes <- c(6, 5, 6, 3)
fig.3b <- ggplot(CszEMM, aes(x = Repro.treatment, y = response, fill = Dev.treatment)) +
geom_jitter(data = ClutchSz,
aes(x = Repro.treatment, y = Clutch.size..1st., colour = Dev.treatment),
position = position_jitterdodge(jitter.width = 0.2, jitter.height = 0, dodge.width = 0.6),
shape = 21, size = 1.5, alpha = 0.2) +
geom_errorbar(aes(ymin = asymp.LCL, ymax = asymp.UCL, 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_sizes, group = interaction(Dev.treatment, Repro.treatment)),
position = position_dodge(0.6), vjust = -1, hjust = -0.75, size = 4, family = "Helvetica", colour = "black") +
labs(x = NULL,
y = "Clutch size") +
scale_y_continuous(breaks = c(150, 250, 350)) +
scale_colour_manual(
name = "Developmental temperature",
values = c("Control" = "blue", "Warm" = "orange")) +
scale_fill_manual(
name = "Developmental temperature",
values = c("Control" = "blue", "Warm" = "orange")) +
theme_classic() +
theme(
text = element_text(family = "Helvetica", colour = "black"),
axis.title = element_text(size = 14, face = "bold", colour = "black"),
axis.text = element_text(size = 12, colour = "black"),
legend.position = "none"
)
fig.3b
Total number of eggs per pair
1. Error structure assessment
poisson.TotEggs <- glmmTMB(Total.number.of.eggs ~ Dev.treatment * Repro.treatment, data=Totals, family="poisson")
negb.TotEggs <- glmmTMB(Total.number.of.eggs ~ Dev.treatment * Repro.treatment, data=Totals, family="nbinom1")res.poisson.TotEggs<-simulateResiduals(poisson.TotEggs, plot=T)
res.negb.TotEggs<-simulateResiduals(negb.TotEggs, plot=T)
Poisson has issues with dispersion. Choosing negative binomial.
2. Maternal effects
Correlation with Female SL
lm.1 = lm(Total.number.of.eggs ~ SL...Female..cm., data=Parental)
summary(lm.1)##
## Call:
## lm(formula = Total.number.of.eggs ~ SL...Female..cm., data = Parental)
##
## Residuals:
## Min 1Q Median 3Q Max
## -486.81 -149.02 30.63 211.06 395.37
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -391.6 1173.8 -0.334 0.743
## SL...Female..cm. 112.6 131.9 0.854 0.405
##
## Residual standard error: 267 on 18 degrees of freedom
## Multiple R-squared: 0.0389, Adjusted R-squared: -0.01449
## F-statistic: 0.7285 on 1 and 18 DF, p-value: 0.4046Correlation with Female W
lm.2 = lm(Total.number.of.eggs ~ Weight...Female..g., data=Parental)
summary(lm.2)##
## Call:
## lm(formula = Total.number.of.eggs ~ Weight...Female..g., data = Parental)
##
## Residuals:
## Min 1Q Median 3Q Max
## -572.3 -135.8 -11.7 212.2 321.3
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.319 377.387 0.014 0.989
## Weight...Female..g. 20.398 12.606 1.618 0.123
##
## Residual standard error: 254.5 on 18 degrees of freedom
## Multiple R-squared: 0.127, Adjusted R-squared: 0.07849
## F-statistic: 2.618 on 1 and 18 DF, p-value: 0.123Correlation with Female FK
lm.3 = lm(Total.number.of.eggs ~ Condition...Female, data=Parental)
summary(lm.3)##
## Call:
## lm(formula = Total.number.of.eggs ~ Condition...Female, data = Parental)
##
## Residuals:
## Min 1Q Median 3Q Max
## -358.83 -269.77 61.99 227.05 340.05
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 70.38 491.94 0.143 0.888
## Condition...Female 127.84 115.92 1.103 0.285
##
## Residual standard error: 263.6 on 18 degrees of freedom
## Multiple R-squared: 0.06329, Adjusted R-squared: 0.01125
## F-statistic: 1.216 on 1 and 18 DF, p-value: 0.2846No maternal traits to be considered for covariate.
3. Final model
Model formula
negb.TotEggs <- glmmTMB(Total.number.of.eggs ~ Dev.treatment * Repro.treatment, data=Totals, family="nbinom1")Performance check
performance::check_zeroinflation(negb.TotEggs)## Model has no observed zeros in the response variable.## NULLperformance::check_model(negb.TotEggs, residual_type = "response")## `check_outliers()` does not yet support models of class `glmmTMB`.
hist(residuals(negb.TotEggs), col="darkgray")
4. Summary statistics … cited in Table S6
summary(negb.TotEggs)## Family: nbinom1 ( log )
## Formula: Total.number.of.eggs ~ Dev.treatment * Repro.treatment
## Data: Totals
##
## AIC BIC logLik deviance df.resid
## 283.0 288.0 -136.5 273.0 15
##
##
## Dispersion parameter for nbinom1 family (): 93.2
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.61711 0.14179 46.67 <2e-16 ***
## Dev.treatmentWarm -0.47756 0.23780 -2.01 0.0446 *
## Repro.treatmentWarm -0.06107 0.20154 -0.30 0.7619
## Dev.treatmentWarm:Repro.treatmentWarm -0.12680 0.38696 -0.33 0.7431
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 15. Anova … cited in Table S3
Anova(negb.TotEggs)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Total.number.of.eggs
## Chisq Df Pr(>Chisq)
## Dev.treatment 7.8089 1 0.005199 **
## Repro.treatment 0.3078 1 0.579016
## Dev.treatment:Repro.treatment 0.1074 1 0.743141
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 16. Emmeans
TEggEMM = (emmeans(negb.TotEggs, ~ Dev.treatment * Repro.treatment, type = "response") %>% as.data.frame)
TEggEMM ## Dev.treatment Repro.treatment response SE df asymp.LCL asymp.UCL
## Control Control 747.7803 106.02737 Inf 566.3470 987.3371
## Warm Control 463.8450 90.03077 Inf 317.0717 678.5600
## Control Warm 703.4829 102.72245 Inf 528.3981 936.5821
## Warm Warm 384.3985 104.39604 Inf 225.7408 654.5660
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale7. Plot … Figure 3c
fig.3c <- ggplot(TEggEMM, aes(x = Repro.treatment, y = response, colour = Dev.treatment)) +
geom_jitter(data = Totals,
aes(x = Repro.treatment, y = Total.number.of.eggs, colour = Dev.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 = asymp.LCL, ymax = asymp.UCL, 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_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 = "Total eggs per pair") +
scale_colour_manual(
name = "Developmental temperature",
values = c("Control" = "blue", "Warm" = "orange")) +
scale_fill_manual(
name = "Developmental temperature",
values = c("Control" = "blue", "Warm" = "orange")) +
theme(
text = element_text(family = "Helvetica", colour = "black"),
axis.title = element_text(size = 14, face = "bold", colour = "black"),
axis.text = element_text(size = 12, colour = "black"),
legend.position = "none"
)
fig.3c
Egg size
1. Maternal effects
Correlation with Female SL
lm.1 = lm(Average.of.Egg.area..mm2. ~ SL...Female..cm., data=Parental)
summary(lm.1)##
## Call:
## lm(formula = Average.of.Egg.area..mm2. ~ SL...Female..cm., data = Parental)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.74644 -0.32307 -0.01939 0.28466 1.17278
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.9405 2.1428 3.239 0.00455 **
## SL...Female..cm. -0.2509 0.2408 -1.042 0.31133
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4875 on 18 degrees of freedom
## Multiple R-squared: 0.05686, Adjusted R-squared: 0.004465
## F-statistic: 1.085 on 1 and 18 DF, p-value: 0.3113Correlation with Female W
lm.2 = lm(Average.of.Egg.area..mm2. ~ Weight...Female..g., data=Parental)
summary(lm.2)##
## Call:
## lm(formula = Average.of.Egg.area..mm2. ~ Weight...Female..g.,
## data = Parental)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.86639 -0.25074 -0.04564 0.27997 1.08511
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.32925 0.73877 5.860 1.5e-05 ***
## Weight...Female..g. 0.01291 0.02468 0.523 0.607
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4982 on 18 degrees of freedom
## Multiple R-squared: 0.01497, Adjusted R-squared: -0.03976
## F-statistic: 0.2735 on 1 and 18 DF, p-value: 0.6074Correlation with Female FK
lm.3 = lm(Average.of.Egg.area..mm2. ~ Condition...Female, data=Parental)
summary(lm.3)##
## Call:
## lm(formula = Average.of.Egg.area..mm2. ~ Condition...Female,
## data = Parental)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.8494 -0.1910 -0.0315 0.2445 0.7583
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.7027 0.8063 3.352 0.00355 **
## Condition...Female 0.4767 0.1900 2.509 0.02189 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4321 on 18 degrees of freedom
## Multiple R-squared: 0.2591, Adjusted R-squared: 0.218
## F-statistic: 6.296 on 1 and 18 DF, p-value: 0.02189Female FK to be considerred as a potential covariate.
2. Covariate testing
lmer.EAr = lmer(Egg.area..mm2. ~ Dev.treatment * Repro.treatment + (1 |Sump/Female.ID), data = EggArea)## boundary (singular) fit: see help('isSingular')lmer.EAra = lmer(Egg.area..mm2. ~ Dev.treatment * Repro.treatment + Female.condition + (1 |Sump/Female.ID), data = EggArea)## boundary (singular) fit: see help('isSingular')AIC(lmer.EAr,lmer.EAra)## df AIC
## lmer.EAr 7 109.6830
## lmer.EAra 8 110.9649Female FK does not improve model fit. Not included in the final
model.
How maternal FK influences thermal effects on egg size is later explored
in the Rmarkdown.
3. Final model
Model formula
lmer.EAr = lmer(Egg.area..mm2. ~ Dev.treatment * Repro.treatment + (1 |Sump/Female.ID), data = EggArea)## boundary (singular) fit: see help('isSingular')Performance check
performance::check_model(lmer.EAr, check="homogeneity") 
performance::check_model(lmer.EAr, check="outliers")
performance::check_model(lmer.EAr, check="qq", detrend = FALSE) 
performance::check_model(lmer.EAr, check="normality") 
performance::check_model(lmer.EAr, check="linearity") 
performance::check_model(lmer.EAr, check="pp_check")
hist(residuals(lmer.EAr), col="darkgray")
outlierTest(lmer.EAr)## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
## rstudent unadjusted p-value Bonferroni p
## 149 2.892581 0.0042645 0.848644. Summary statistics … cited in Table S6
summary(lmer.EAr)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## Egg.area..mm2. ~ Dev.treatment * Repro.treatment + (1 | Sump/Female.ID)
## Data: EggArea
##
## REML criterion at convergence: 95.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.59812 -0.58923 -0.04186 0.61334 2.75039
##
## Random effects:
## Groups Name Variance Std.Dev.
## Female.ID:Sump (Intercept) 0.13324 0.3650
## Sump (Intercept) 0.00000 0.0000
## Residual 0.06903 0.2627
## Number of obs: 199, groups: Female.ID:Sump, 20; Sump, 6
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 5.1105 0.1528 15.9992 33.439
## Dev.treatmentWarm -0.2587 0.2267 15.9992 -1.141
## Repro.treatmentWarm -0.8411 0.2162 16.0136 -3.891
## Dev.treatmentWarm:Repro.treatmentWarm 0.5513 0.3485 16.0047 1.582
## Pr(>|t|)
## (Intercept) 3.1e-16 ***
## Dev.treatmentWarm 0.2706
## Repro.treatmentWarm 0.0013 **
## Dev.treatmentWarm:Repro.treatmentWarm 0.1332
## ---
## 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.707 0.477
## Dv.trtW:R.W 0.438 -0.650 -0.620
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')5. anova … cited in Table S3
anova(lmer.EAr)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Dev.treatment 0.00066 0.00066 1 16.005 0.0095 0.923416
## Repro.treatment 0.72674 0.72674 1 16.005 10.5278 0.005076 **
## Dev.treatment:Repro.treatment 0.17274 0.17274 1 16.005 2.5023 0.133235
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 16. emmeans
EaEMM = (emmeans(lmer.EAr, ~ Dev.treatment * Repro.treatment, type = "response") %>% as.data.frame)
EaEMM## Dev.treatment Repro.treatment emmean SE df lower.CL upper.CL
## Control Control 5.110450 0.1920384 3.79 4.565245 5.655655
## Warm Control 4.851800 0.1715306 9.21 4.465087 5.238513
## Control Warm 4.269346 0.1641498 6.10 3.869275 4.669417
## Warm Warm 4.562033 0.2444243 10.19 4.018796 5.105270
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.957. Plot … Figure 3d
sample_size_eggs <- c(60, 50, 59, 30)
fig.3d <- ggplot(EaEMM, aes(x = Repro.treatment, y = emmean, colour = Dev.treatment)) +
geom_jitter(data = EggArea,
aes(x = Repro.treatment, y = Egg.area..mm2., colour = Dev.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_eggs, 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 = "Egg size (mm²)") +
scale_y_continuous(breaks = c(3.5, 4.5, 5.5, 6.5)) +
scale_colour_manual(
name = "Developmental temperature",
values = c("Control" = "blue", "Warm" = "orange")) +
scale_fill_manual(
name = "Developmental temperature",
values = c("Control" = "blue", "Warm" = "orange")) +
theme(
text = element_text(family = "Helvetica"),
axis.text = element_text(size = 12, colour = "black"),
axis.title = element_text(size = 14, face = "bold", colour = "black"),
legend.position = "none"
)
fig.3d
How matenral covariates impacted thermal effects
For traits significantly correlated with maternal traits (clutch size and Female SL, egg size and Female FK), the summary statistics and anova are compared between models with and without the selected maternal traits to evaluate their impact on observed thermal effects.
Clutch size and maternal SL
Summary statistics
Summary statistics for the model without Female SL are provided
here.
Please refer to the final model for summary statistics of models
including Female SL.
negb.Csz <- glmmTMB(Clutch.size..1st. ~ Dev.treatment * Repro.treatment, data=ClutchSz, family="nbinom1")
summary(negb.Csz)## Family: nbinom1 ( log )
## Formula: Clutch.size..1st. ~ Dev.treatment * Repro.treatment
## Data: ClutchSz
##
## AIC BIC logLik deviance df.resid
## 231.3 236.3 -110.7 221.3 15
##
##
## Dispersion parameter for nbinom1 family (): 14.1
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.56238 0.09739 57.11 <2e-16 ***
## Dev.treatmentWarm -0.15325 0.15003 -1.02 0.307
## Repro.treatmentWarm 0.11465 0.13352 0.86 0.391
## Dev.treatmentWarm:Repro.treatmentWarm -0.03116 0.22588 -0.14 0.890
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Anova
Final model (Maternal SL added):
Anova(negb.Csza)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Clutch.size..1st.
## Chisq Df Pr(>Chisq)
## Dev.treatment 0.5744 1 0.44852
## Repro.treatment 0.2607 1 0.60965
## Female.SL..cm. 6.0207 1 0.01414 *
## Dev.treatment:Repro.treatment 0.2962 1 0.58630
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Model without Maternal SL
Anova(negb.Csz)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Clutch.size..1st.
## Chisq Df Pr(>Chisq)
## Dev.treatment 2.2162 1 0.1366
## Repro.treatment 0.9283 1 0.3353
## Dev.treatment:Repro.treatment 0.0190 1 0.8903When maternal SL is not considered, developmental warming generally
reduces clutch size. However, this effect diminishes in the final model
when maternal SL is included, suggesting the impact of developmental
warming is manifested through its impact on adult size.
Adults exposed to warming during development were smaller, leading to
smaller clutches.
Egg size and Female FK
Summary statistics
Summary statistics for the model with Female FK are provided
here.
Please refer above for the final model for summary statistics of models
without Female FK.
lmer.EAra = lmer(Egg.area..mm2. ~ Dev.treatment * Repro.treatment + Female.condition + (1 |Sump/Female.ID), data = EggArea)## boundary (singular) fit: see help('isSingular')summary(lmer.EAra)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Egg.area..mm2. ~ Dev.treatment * Repro.treatment + Female.condition +
## (1 | Sump/Female.ID)
## Data: EggArea
##
## REML criterion at convergence: 95
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.61436 -0.58453 -0.05911 0.61648 2.77670
##
## Random effects:
## Groups Name Variance Std.Dev.
## Female.ID:Sump (Intercept) 0.12208 0.3494
## Sump (Intercept) 0.00000 0.0000
## Residual 0.06903 0.2627
## Number of obs: 199, groups: Female.ID:Sump, 20; Sump, 6
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 3.8544 0.8267 14.9992 4.662
## Dev.treatmentWarm -0.1501 0.2286 14.9991 -0.657
## Repro.treatmentWarm -0.7011 0.2264 15.0120 -3.098
## Female.condition 0.2756 0.1785 14.9992 1.544
## Dev.treatmentWarm:Repro.treatmentWarm 0.4756 0.3380 15.0048 1.407
## Pr(>|t|)
## (Intercept) 0.000307 ***
## Dev.treatmentWarm 0.521363
## Repro.treatmentWarm 0.007347 **
## Female.condition 0.143480
## Dev.treatmentWarm:Repro.treatmentWarm 0.179765
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Dv.trW Rpr.tW Fml.cn
## Dv.trtmntWr -0.417
## Rpr.trtmntW -0.509 0.539
## Femal.cndtn -0.984 0.308 0.401
## Dv.trtW:R.W 0.220 -0.657 -0.620 -0.145
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')Singular fit as no variance is assigned to Sump.
Anova
Final model (without Female FK)
anova(lmer.EAr)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Dev.treatment 0.00066 0.00066 1 16.005 0.0095 0.923416
## Repro.treatment 0.72674 0.72674 1 16.005 10.5278 0.005076 **
## Dev.treatment:Repro.treatment 0.17274 0.17274 1 16.005 2.5023 0.133235
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Model with Female FK
lmer.EAra = lmer(Egg.area..mm2. ~ Dev.treatment * Repro.treatment + Female.condition + (1 |Sump/Female.ID), data = EggArea)## boundary (singular) fit: see help('isSingular')anova(lmer.EAra)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Dev.treatment 0.01767 0.01767 1 15.004 0.2559 0.62028
## Repro.treatment 0.45854 0.45854 1 15.004 6.6426 0.02102 *
## Female.condition 0.16451 0.16451 1 14.999 2.3832 0.14348
## Dev.treatment:Repro.treatment 0.13668 0.13668 1 15.005 1.9800 0.17976
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Female FK does not change the overall thermal impact.
Combining plots for Figure 3
fig.3_legend <- fig.3a + theme(legend.position = "top")
legend <- get_legend(fig.3_legend)
fig.Repro <- ggarrange(
fig.3a, fig.3b, fig.3c, fig.3d,
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.Repro_with_legend <- ggarrange(
legend, fig.Repro,
ncol = 1, heights = c(0.1, 0.9)
)
fig.Repro_with_legend <- annotate_figure(
fig.Repro_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 3.pdf",
plot = fig.Repro_with_legend,
width = 169,
height = 180,
units = "mm",
dpi = 600
)