Chapter 1 _ MHW impacts on juveniles
J Donelson and Y Yasutake
Created Jan 2024, Updated Jan 2026
R Markdown
library(ggplot2)
library(tidyverse)
library(lme4)
library(lmerTest)
library(emmeans)
library(car)
library(performance)
library(readxl)Data loading
Egg <- read_excel(
path = "Embryonic exposure to a marine heatwave restricts juvenile growth in a coral reef damselfish.xlsx",
sheet = "Egg size",
trim_ws = TRUE)
Egg$Treatment = factor(Egg$Treatment)
Egg$Clutch = factor(Egg$Clutch)
Egg$Parental.tank = factor(Egg$Parental.tank)
Egg$Male = factor(Egg$Male)
Egg$Female = factor(Egg$Female)
Egg$Sample = factor(Egg$Sample)
Egg$Egg.size = as.numeric(Egg$Egg.size)
str(Egg)## tibble [170 × 13] (S3: tbl_df/tbl/data.frame)
## $ Clutch : Factor w/ 17 levels "61","92","94",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Treatment : Factor w/ 2 levels "Control","Embryo": 1 1 1 1 1 1 1 1 1 1 ...
## $ Parental.tank: Factor w/ 15 levels "14","17","30",..: 14 14 14 14 14 14 14 14 14 14 ...
## $ Male : Factor w/ 15 levels "CARL217","CARL335",..: 15 15 15 15 15 15 15 15 15 15 ...
## $ Female : Factor w/ 15 levels "CARL226","CARL345",..: 14 14 14 14 14 14 14 14 14 14 ...
## $ MaleSL : num [1:170] 9.71 9.71 9.71 9.71 9.71 ...
## $ MaleW : num [1:170] 46 46 46 46 46 ...
## $ MaleFK : num [1:170] 5.02 5.02 5.02 5.02 5.02 ...
## $ FemaleSL : num [1:170] 9.66 9.66 9.66 9.66 9.66 ...
## $ FemaleW : num [1:170] 43.9 43.9 43.9 43.9 43.9 ...
## $ FemaleFK : num [1:170] 4.88 4.88 4.88 4.88 4.88 ...
## $ Sample : Factor w/ 10 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ Egg.size : num [1:170] 4.9 5.1 4.8 5.6 5.2 5 5.3 5.1 5.3 5 ...MHW <- read_excel(
path = "Embryonic exposure to a marine heatwave restricts juvenile growth in a coral reef damselfish.xlsx",
sheet = "Day 0-60",
trim_ws = TRUE)
MHW$Clutch = factor(MHW$Clutch)
MHW$Treatment = factor(MHW$Treatment)
MHW$Day = factor(MHW$Day)
MHW$Male = factor(MHW$Male)
MHW$Female = factor(MHW$Female)
MHW$Fish.rep = factor(MHW$Fish.rep)
MHW$Parental.tank = factor(MHW$Parental.tank)
MHW$Tank = factor(MHW$Tank)
MHW$SL = as.numeric(MHW$SL)
MHW$W = as.numeric(MHW$W)
MHW$YA = as.numeric(MHW$YA)
MHW$FemaleSL = as.numeric(MHW$FemaleSL)
MHW$FemaleW = as.numeric(MHW$FemaleW)
MHW$FemaleFK = as.numeric(MHW$FemaleFK)
MHW$MaleSL = as.numeric(MHW$MaleSL)
MHW$MaleW = as.numeric(MHW$MaleW)
MHW$MaleFK = as.numeric(MHW$MaleFK)
MHW$LogFemaleSL <- log(MHW$FemaleSL)
MHW$LogMaleSL <- log(MHW$MaleSL)
MHW$LogFemaleFK <- log(MHW$FemaleFK)
MHW$LogMaleFK <- log(MHW$MaleFK)
MHW$LogSL <- log(MHW$SL)
MHW$LogW <- log(MHW$W)
str(MHW)## tibble [2,304 × 24] (S3: tbl_df/tbl/data.frame)
## $ Clutch : Factor w/ 17 levels "61","92","94",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Treatment : Factor w/ 4 levels "Control","Early exposure",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Parental.tank: Factor w/ 15 levels "14","17","30",..: 14 14 14 14 14 14 14 14 14 14 ...
## $ Male : Factor w/ 15 levels "CARL217","CARL335",..: 15 15 15 15 15 15 15 15 15 15 ...
## $ Female : Factor w/ 15 levels "CARL226","CARL345",..: 14 14 14 14 14 14 14 14 14 14 ...
## $ MaleSL : num [1:2304] 9.71 9.71 9.71 9.71 9.71 ...
## $ MaleW : num [1:2304] 46 46 46 46 46 ...
## $ MaleFK : num [1:2304] 5.02 5.02 5.02 5.02 5.02 ...
## $ FemaleSL : num [1:2304] 9.66 9.66 9.66 9.66 9.66 ...
## $ FemaleW : num [1:2304] 43.9 43.9 43.9 43.9 43.9 ...
## $ FemaleFK : num [1:2304] 4.88 4.88 4.88 4.88 4.88 ...
## $ Day : Factor w/ 5 levels "0","15","30",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Tank : Factor w/ 76 levels "14","17","30",..: 59 59 59 59 59 59 59 59 59 59 ...
## $ Fish.rep : Factor w/ 30 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ SL : num [1:2304] 0.513 0.518 0.532 0.534 0.535 ...
## $ W : num [1:2304] 0.0028 0.0026 0.0024 0.0029 0.0029 0.0028 0.0022 0.0029 0.0029 0.0028 ...
## $ YA : num [1:2304] 0.016 0.012 0.014 0.015 0.013 0.011 0.011 0.014 0.014 0.013 ...
## $ Note : chr [1:2304] NA NA NA NA ...
## $ LogFemaleSL : num [1:2304] 2.27 2.27 2.27 2.27 2.27 ...
## $ LogMaleSL : num [1:2304] 2.27 2.27 2.27 2.27 2.27 ...
## $ LogFemaleFK : num [1:2304] 1.58 1.58 1.58 1.58 1.58 ...
## $ LogMaleFK : num [1:2304] 1.61 1.61 1.61 1.61 1.61 ...
## $ LogSL : num [1:2304] -0.667 -0.657 -0.631 -0.627 -0.626 ...
## $ LogW : num [1:2304] -5.88 -5.95 -6.03 -5.84 -5.84 ...MHWembryo <- MHW %>%
filter(Treatment != "Early exposure",
Treatment != "Late exposure",
Day != 0)
str(MHWembryo)## tibble [1,182 × 24] (S3: tbl_df/tbl/data.frame)
## $ Clutch : Factor w/ 17 levels "61","92","94",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Treatment : Factor w/ 4 levels "Control","Early exposure",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Parental.tank: Factor w/ 15 levels "14","17","30",..: 14 14 14 14 14 14 14 14 14 14 ...
## $ Male : Factor w/ 15 levels "CARL217","CARL335",..: 15 15 15 15 15 15 15 15 15 15 ...
## $ Female : Factor w/ 15 levels "CARL226","CARL345",..: 14 14 14 14 14 14 14 14 14 14 ...
## $ MaleSL : num [1:1182] 9.71 9.71 9.71 9.71 9.71 ...
## $ MaleW : num [1:1182] 46 46 46 46 46 ...
## $ MaleFK : num [1:1182] 5.02 5.02 5.02 5.02 5.02 ...
## $ FemaleSL : num [1:1182] 9.66 9.66 9.66 9.66 9.66 ...
## $ FemaleW : num [1:1182] 43.9 43.9 43.9 43.9 43.9 ...
## $ FemaleFK : num [1:1182] 4.88 4.88 4.88 4.88 4.88 ...
## $ Day : Factor w/ 5 levels "0","15","30",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ Tank : Factor w/ 76 levels "14","17","30",..: 16 16 16 16 16 16 16 16 16 16 ...
## $ Fish.rep : Factor w/ 30 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ SL : num [1:1182] 0.892 0.949 0.89 0.861 0.902 ...
## $ W : num [1:1182] 0.0283 0.0355 0.0251 0.0244 0.0266 0.0228 0.0292 0.0268 0.0316 0.0306 ...
## $ YA : num [1:1182] NA NA NA NA NA NA NA NA NA NA ...
## $ Note : chr [1:1182] NA NA NA NA ...
## $ LogFemaleSL : num [1:1182] 2.27 2.27 2.27 2.27 2.27 ...
## $ LogMaleSL : num [1:1182] 2.27 2.27 2.27 2.27 2.27 ...
## $ LogFemaleFK : num [1:1182] 1.58 1.58 1.58 1.58 1.58 ...
## $ LogMaleFK : num [1:1182] 1.61 1.61 1.61 1.61 1.61 ...
## $ LogSL : num [1:1182] -0.114 -0.0524 -0.1162 -0.1496 -0.103 ...
## $ LogW : num [1:1182] -3.56 -3.34 -3.68 -3.71 -3.63 ...MHWday15_30 <- MHW %>%
filter(Day %in% c("15", "30"),
Treatment != "Embryo",
!Parental.tank %in% c("31", "144", "216")) # removing pairs that are not present in all treatments
str(MHWday15_30)## tibble [444 × 24] (S3: tbl_df/tbl/data.frame)
## $ Clutch : Factor w/ 17 levels "61","92","94",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Treatment : Factor w/ 4 levels "Control","Early exposure",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Parental.tank: Factor w/ 15 levels "14","17","30",..: 14 14 14 14 14 14 14 14 14 14 ...
## $ Male : Factor w/ 15 levels "CARL217","CARL335",..: 15 15 15 15 15 15 15 15 15 15 ...
## $ Female : Factor w/ 15 levels "CARL226","CARL345",..: 14 14 14 14 14 14 14 14 14 14 ...
## $ MaleSL : num [1:444] 9.71 9.71 9.71 9.71 9.71 ...
## $ MaleW : num [1:444] 46 46 46 46 46 ...
## $ MaleFK : num [1:444] 5.02 5.02 5.02 5.02 5.02 ...
## $ FemaleSL : num [1:444] 9.66 9.66 9.66 9.66 9.66 ...
## $ FemaleW : num [1:444] 43.9 43.9 43.9 43.9 43.9 ...
## $ FemaleFK : num [1:444] 4.88 4.88 4.88 4.88 4.88 ...
## $ Day : Factor w/ 5 levels "0","15","30",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ Tank : Factor w/ 76 levels "14","17","30",..: 16 16 16 16 16 16 16 16 16 16 ...
## $ Fish.rep : Factor w/ 30 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ SL : num [1:444] 0.892 0.949 0.89 0.861 0.902 ...
## $ W : num [1:444] 0.0283 0.0355 0.0251 0.0244 0.0266 0.0228 0.0292 0.0268 0.0316 0.0306 ...
## $ YA : num [1:444] NA NA NA NA NA NA NA NA NA NA ...
## $ Note : chr [1:444] NA NA NA NA ...
## $ LogFemaleSL : num [1:444] 2.27 2.27 2.27 2.27 2.27 ...
## $ LogMaleSL : num [1:444] 2.27 2.27 2.27 2.27 2.27 ...
## $ LogFemaleFK : num [1:444] 1.58 1.58 1.58 1.58 1.58 ...
## $ LogMaleFK : num [1:444] 1.61 1.61 1.61 1.61 1.61 ...
## $ LogSL : num [1:444] -0.114 -0.0524 -0.1162 -0.1496 -0.103 ...
## $ LogW : num [1:444] -3.56 -3.34 -3.68 -3.71 -3.63 ...MHWday45_60 <- MHW %>%
filter(Day %in% c("45", "60"),
Treatment != "Embryo",
!Parental.tank %in% c("31", "144", "216")) # removing pairs that are not present in all treatments
str(MHWday45_60)## tibble [602 × 24] (S3: tbl_df/tbl/data.frame)
## $ Clutch : Factor w/ 17 levels "61","92","94",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Treatment : Factor w/ 4 levels "Control","Early exposure",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Parental.tank: Factor w/ 15 levels "14","17","30",..: 14 14 14 14 14 14 14 14 14 14 ...
## $ Male : Factor w/ 15 levels "CARL217","CARL335",..: 15 15 15 15 15 15 15 15 15 15 ...
## $ Female : Factor w/ 15 levels "CARL226","CARL345",..: 14 14 14 14 14 14 14 14 14 14 ...
## $ MaleSL : num [1:602] 9.71 9.71 9.71 9.71 9.71 ...
## $ MaleW : num [1:602] 46 46 46 46 46 ...
## $ MaleFK : num [1:602] 5.02 5.02 5.02 5.02 5.02 ...
## $ FemaleSL : num [1:602] 9.66 9.66 9.66 9.66 9.66 ...
## $ FemaleW : num [1:602] 43.9 43.9 43.9 43.9 43.9 ...
## $ FemaleFK : num [1:602] 4.88 4.88 4.88 4.88 4.88 ...
## $ Day : Factor w/ 5 levels "0","15","30",..: 4 4 4 4 4 4 4 4 4 4 ...
## $ Tank : Factor w/ 76 levels "14","17","30",..: 39 39 39 39 39 39 39 39 39 39 ...
## $ Fish.rep : Factor w/ 30 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ SL : num [1:602] 2.03 1.73 1.48 1.81 1.94 ...
## $ W : num [1:602] 0.314 0.209 0.111 0.242 0.276 ...
## $ YA : num [1:602] NA NA NA NA NA NA NA NA NA NA ...
## $ Note : chr [1:602] NA NA NA NA ...
## $ LogFemaleSL : num [1:602] 2.27 2.27 2.27 2.27 2.27 ...
## $ LogMaleSL : num [1:602] 2.27 2.27 2.27 2.27 2.27 ...
## $ LogFemaleFK : num [1:602] 1.58 1.58 1.58 1.58 1.58 ...
## $ LogMaleFK : num [1:602] 1.61 1.61 1.61 1.61 1.61 ...
## $ LogSL : num [1:602] 0.707 0.545 0.39 0.593 0.662 ...
## $ LogW : num [1:602] -1.16 -1.57 -2.2 -1.42 -1.29 ...Hatching = MHW %>%
filter(Day %in% c("0"))
str(Hatching)## tibble [469 × 24] (S3: tbl_df/tbl/data.frame)
## $ Clutch : Factor w/ 17 levels "61","92","94",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Treatment : Factor w/ 4 levels "Control","Early exposure",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Parental.tank: Factor w/ 15 levels "14","17","30",..: 14 14 14 14 14 14 14 14 14 14 ...
## $ Male : Factor w/ 15 levels "CARL217","CARL335",..: 15 15 15 15 15 15 15 15 15 15 ...
## $ Female : Factor w/ 15 levels "CARL226","CARL345",..: 14 14 14 14 14 14 14 14 14 14 ...
## $ MaleSL : num [1:469] 9.71 9.71 9.71 9.71 9.71 ...
## $ MaleW : num [1:469] 46 46 46 46 46 ...
## $ MaleFK : num [1:469] 5.02 5.02 5.02 5.02 5.02 ...
## $ FemaleSL : num [1:469] 9.66 9.66 9.66 9.66 9.66 ...
## $ FemaleW : num [1:469] 43.9 43.9 43.9 43.9 43.9 ...
## $ FemaleFK : num [1:469] 4.88 4.88 4.88 4.88 4.88 ...
## $ Day : Factor w/ 5 levels "0","15","30",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Tank : Factor w/ 76 levels "14","17","30",..: 59 59 59 59 59 59 59 59 59 59 ...
## $ Fish.rep : Factor w/ 30 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ SL : num [1:469] 0.513 0.518 0.532 0.534 0.535 ...
## $ W : num [1:469] 0.0028 0.0026 0.0024 0.0029 0.0029 0.0028 0.0022 0.0029 0.0029 0.0028 ...
## $ YA : num [1:469] 0.016 0.012 0.014 0.015 0.013 0.011 0.011 0.014 0.014 0.013 ...
## $ Note : chr [1:469] NA NA NA NA ...
## $ LogFemaleSL : num [1:469] 2.27 2.27 2.27 2.27 2.27 ...
## $ LogMaleSL : num [1:469] 2.27 2.27 2.27 2.27 2.27 ...
## $ LogFemaleFK : num [1:469] 1.58 1.58 1.58 1.58 1.58 ...
## $ LogMaleFK : num [1:469] 1.61 1.61 1.61 1.61 1.61 ...
## $ LogSL : num [1:469] -0.667 -0.657 -0.631 -0.627 -0.626 ...
## $ LogW : num [1:469] -5.88 -5.95 -6.03 -5.84 -5.84 ...EarlyExposure <- MHW %>%
filter(Day %in% c("15", "30", "45", "60"),
Treatment %in% c("Control", "Early exposure"),
!Parental.tank %in% c("31", "144", "216"))
str(EarlyExposure)## tibble [825 × 24] (S3: tbl_df/tbl/data.frame)
## $ Clutch : Factor w/ 17 levels "61","92","94",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Treatment : Factor w/ 4 levels "Control","Early exposure",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Parental.tank: Factor w/ 15 levels "14","17","30",..: 14 14 14 14 14 14 14 14 14 14 ...
## $ Male : Factor w/ 15 levels "CARL217","CARL335",..: 15 15 15 15 15 15 15 15 15 15 ...
## $ Female : Factor w/ 15 levels "CARL226","CARL345",..: 14 14 14 14 14 14 14 14 14 14 ...
## $ MaleSL : num [1:825] 9.71 9.71 9.71 9.71 9.71 ...
## $ MaleW : num [1:825] 46 46 46 46 46 ...
## $ MaleFK : num [1:825] 5.02 5.02 5.02 5.02 5.02 ...
## $ FemaleSL : num [1:825] 9.66 9.66 9.66 9.66 9.66 ...
## $ FemaleW : num [1:825] 43.9 43.9 43.9 43.9 43.9 ...
## $ FemaleFK : num [1:825] 4.88 4.88 4.88 4.88 4.88 ...
## $ Day : Factor w/ 5 levels "0","15","30",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ Tank : Factor w/ 76 levels "14","17","30",..: 16 16 16 16 16 16 16 16 16 16 ...
## $ Fish.rep : Factor w/ 30 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ SL : num [1:825] 0.892 0.949 0.89 0.861 0.902 ...
## $ W : num [1:825] 0.0283 0.0355 0.0251 0.0244 0.0266 0.0228 0.0292 0.0268 0.0316 0.0306 ...
## $ YA : num [1:825] NA NA NA NA NA NA NA NA NA NA ...
## $ Note : chr [1:825] NA NA NA NA ...
## $ LogFemaleSL : num [1:825] 2.27 2.27 2.27 2.27 2.27 ...
## $ LogMaleSL : num [1:825] 2.27 2.27 2.27 2.27 2.27 ...
## $ LogFemaleFK : num [1:825] 1.58 1.58 1.58 1.58 1.58 ...
## $ LogMaleFK : num [1:825] 1.61 1.61 1.61 1.61 1.61 ...
## $ LogSL : num [1:825] -0.114 -0.0524 -0.1162 -0.1496 -0.103 ...
## $ LogW : num [1:825] -3.56 -3.34 -3.68 -3.71 -3.63 ...Overview
For each trait, information are given as follows:
1. Random effect structures
2. Covariate exploration
3. Final model
4. Model performance
5. Model outputs (summary statistics, anova, emmeans, and pairwise
comparisons where appropriate)
6. Figure
Effects of MHWs during embryonic development
Egg size (Control vs Embryo)
1. Random effect structure
lmer.egg = lmer(Egg.size ~ Treatment + (1|Parental.tank/Clutch), data = Egg)Clutch nested within Parental tank as 2 clutches from Pair 32 and 66 are used (1 each in Control and Embryo).
2. Covariate exploration
lm = lm(Egg.size ~ FemaleSL, data=Egg)
summary(lm)##
## Call:
## lm(formula = Egg.size ~ FemaleSL, data = Egg)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.62273 -0.31047 -0.00942 0.39023 1.50977
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.41824 0.47033 3.015 0.00296 **
## FemaleSL 0.35903 0.04996 7.186 2.08e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.576 on 168 degrees of freedom
## Multiple R-squared: 0.2351, Adjusted R-squared: 0.2305
## F-statistic: 51.63 on 1 and 168 DF, p-value: 2.08e-11lm = lm(Egg.size ~ FemaleFK, data=Egg)
summary(lm)##
## Call:
## lm(formula = Egg.size ~ FemaleFK, data = Egg)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.6907 -0.2640 0.0966 0.3881 1.2655
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.5616 0.5608 1.001 0.318
## FemaleFK 0.8514 0.1128 7.550 2.64e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5691 on 168 degrees of freedom
## Multiple R-squared: 0.2533, Adjusted R-squared: 0.2489
## F-statistic: 57 on 1 and 168 DF, p-value: 2.639e-12Stronger correlation with maternal FK. Testing model fit.
lmer.egg = lmer(Egg.size ~ Treatment + (1|Parental.tank/Clutch), data = Egg)
lmer.egga = lmer(Egg.size ~ Treatment + FemaleFK + (1|Parental.tank/Clutch), data = Egg)
AIC(lmer.egg, lmer.egga)## df AIC
## lmer.egg 5 193.0985
## lmer.egga 6 190.5055Choosing a (with FemaleFK)
3. Final model
lmer.egg = lmer(Egg.size ~ Treatment + FemaleFK + (1|Parental.tank/Clutch), data = Egg)4. Performance check
performance::check_model(lmer.egg, check="homogeneity") 
performance::check_model(lmer.egg, check="outliers") 
performance::check_model(lmer.egg, check="qq", detrend = FALSE) 
performance::check_model(lmer.egg, check="normality") 
performance::check_model(lmer.egg, check="linearity") 
performance::check_model(lmer.egg, check="pp_check")
hist(residuals(lmer.egg), col="darkgray")
shapiro.test(residuals(lmer.egg))##
## Shapiro-Wilk normality test
##
## data: residuals(lmer.egg)
## W = 0.99427, p-value = 0.751qqline(resid(lmer.egg))
outlierTest(lmer.egg)## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
## rstudent unadjusted p-value Bonferroni p
## 92 -3.173765 0.0017991 0.30584#92 identified as outlier but the egg is normal (e.g., not squished). Kept.
5. Outputs
Summary statistics
summary(lmer.egg)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Egg.size ~ Treatment + FemaleFK + (1 | Parental.tank/Clutch)
## Data: Egg
##
## REML criterion at convergence: 178.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.01829 -0.58358 -0.03524 0.55607 2.49202
##
## Random effects:
## Groups Name Variance Std.Dev.
## Clutch:Parental.tank (Intercept) 0.0493 0.2220
## Parental.tank (Intercept) 0.1922 0.4384
## Residual 0.1241 0.3523
## Number of obs: 170, groups: Clutch:Parental.tank, 17; Parental.tank, 15
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.9345 1.6577 13.3030 0.564 0.5823
## TreatmentEmbryo -0.1127 0.1888 3.6123 -0.597 0.5858
## FemaleFK 0.7808 0.3310 13.0192 2.359 0.0346 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) TrtmnE
## TrtmntEmbry -0.225
## FemaleFK -0.996 0.175anova
Anova(lmer.egg, type = 3)## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: Egg.size
## Chisq Df Pr(>Chisq)
## (Intercept) 0.3178 1 0.57293
## Treatment 0.3567 1 0.55035
## FemaleFK 5.5630 1 0.01834 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Emmeans
EmmEgg = (emmeans(lmer.egg, ~ Treatment) %>% as.data.frame)
EmmEgg## Treatment emmean SE df lower.CL upper.CL
## Control 4.805578 0.1758959 12.85 4.425126 5.186030
## Embryo 4.692833 0.1909206 11.83 4.276177 5.109489
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.956. Figure
fig.Egg <- ggplot(EmmEgg,
aes(x = Treatment, y = emmean, fill = Treatment)) +
geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL),
position = position_dodge(width = 0.6),
width = 0.1, colour = "black", linewidth = 0.3) +
geom_point(shape = 22, size = 8, colour = "black", stroke = 0.3,
position = position_dodge(width = 0.6)) +
theme_classic() +
labs(y = "Egg size (mm²)") +
scale_fill_manual(values = c("Control" = "#0072B2",
"Embryo" = "#B2182B")) +
scale_y_continuous(
limits = c(4, 5.5),
breaks = c(4, 4.5, 5, 5.5),
labels = c(4, 4.5, 5, 5.5)) +
theme(
text = element_text(family = "Helvetica"),
axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.title.y = element_text(size = 14, face = "bold", colour = "black"),
axis.text.y = element_text(size = 12, colour = "black"),
legend.position = "none"
)
fig.Egg
Juvenile morphometrics at hatching
Standard length at hatching
1. Random effect structures
lmer.HatchSL = lmer(SL ~ Treatment + (1|Clutch), data = Hatching)
lmer.HatchSLa = lmer(SL ~ Treatment + (1|Parental.tank/Clutch), data = Hatching)
AIC(lmer.HatchSL, lmer.HatchSLa)## df AIC
## lmer.HatchSL 4 -2185.504
## lmer.HatchSLa 5 -2183.892Choosing “a” to account for cases where clutches from the same pair is used.
2. Covariate exploration
lm = lm(SL ~ FemaleSL, data=Hatching)
summary(lm)##
## Call:
## lm(formula = SL ~ FemaleSL, data = Hatching)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.100740 -0.016650 -0.000944 0.016981 0.087383
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.431337 0.014282 30.20 < 2e-16 ***
## FemaleSL 0.011868 0.001508 7.87 2.56e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02885 on 461 degrees of freedom
## (6 observations deleted due to missingness)
## Multiple R-squared: 0.1184, Adjusted R-squared: 0.1165
## F-statistic: 61.93 on 1 and 461 DF, p-value: 2.556e-14lm = lm(SL ~ FemaleFK, data=Hatching)
summary(lm)##
## Call:
## lm(formula = SL ~ FemaleFK, data = Hatching)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.096227 -0.016336 0.000423 0.016889 0.090269
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.496791 0.018663 26.619 <2e-16 ***
## FemaleFK 0.009311 0.003731 2.496 0.0129 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03052 on 461 degrees of freedom
## (6 observations deleted due to missingness)
## Multiple R-squared: 0.01333, Adjusted R-squared: 0.01119
## F-statistic: 6.229 on 1 and 461 DF, p-value: 0.01292lm = lm(SL ~ MaleSL, data=Hatching)
summary(lm)##
## Call:
## lm(formula = SL ~ MaleSL, data = Hatching)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.10040 -0.01730 -0.00104 0.01799 0.09384
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.463214 0.019622 23.607 < 2e-16 ***
## MaleSL 0.008440 0.002064 4.089 5.12e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03019 on 461 degrees of freedom
## (6 observations deleted due to missingness)
## Multiple R-squared: 0.03499, Adjusted R-squared: 0.0329
## F-statistic: 16.72 on 1 and 461 DF, p-value: 5.121e-05lm = lm(SL ~ MaleFK, data=Hatching)
summary(lm)##
## Call:
## lm(formula = SL ~ MaleFK, data = Hatching)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.09737 -0.01724 -0.00031 0.01926 0.08685
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.439137 0.018224 24.097 < 2e-16 ***
## MaleFK 0.020654 0.003605 5.729 1.83e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02969 on 461 degrees of freedom
## (6 observations deleted due to missingness)
## Multiple R-squared: 0.06646, Adjusted R-squared: 0.06443
## F-statistic: 32.82 on 1 and 461 DF, p-value: 1.829e-08Strongest correlation with Female SL. Testing model fit.
lmer.HatchSLa = lmer(SL ~ Treatment + (1|Parental.tank/Clutch), data = Hatching)
lmer.HatchSLb = lmer(SL ~ Treatment + FemaleSL + (1|Parental.tank/Clutch), data = Hatching)
AIC(lmer.HatchSLa, lmer.HatchSLb)## df AIC
## lmer.HatchSLa 5 -2183.892
## lmer.HatchSLb 6 -2176.945Choosing a (without covariate).
3. Final model
lmer.HatchSL = lmer(SL ~ Treatment + (1|Parental.tank/Clutch), data = Hatching)4. Model performance
performance::check_model(lmer.HatchSL, check="homogeneity") 
performance::check_model(lmer.HatchSL, check="outliers") 
performance::check_model(lmer.HatchSL, check="qq", detrend = FALSE) 
performance::check_model(lmer.HatchSL, check="normality") 
performance::check_model(lmer.HatchSL, check="linearity") 
performance::check_model(lmer.HatchSL, check="pp_check", re_formula=NA)
hist(residuals(lmer.HatchSL), col="darkgray")
shapiro.test(residuals(lmer.HatchSL))##
## Shapiro-Wilk normality test
##
## data: residuals(lmer.HatchSL)
## W = 0.98168, p-value = 1.388e-05qqline(resid(lmer.HatchSL))
outlierTest(lmer.HatchSL)## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
## rstudent unadjusted p-value Bonferroni p
## 320 -3.465422 0.00057937 0.26825#320 identified as an outlier but retained as it appears within the expected range for that clutch.
5. Output
Summary statistics
summary(lmer.HatchSL)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: SL ~ Treatment + (1 | Parental.tank/Clutch)
## Data: Hatching
##
## REML criterion at convergence: -2193.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.4067 -0.5318 0.0392 0.6867 2.2585
##
## Random effects:
## Groups Name Variance Std.Dev.
## Clutch:Parental.tank (Intercept) 0.0003737 0.01933
## Parental.tank (Intercept) 0.0001745 0.01321
## Residual 0.0004395 0.02097
## Number of obs: 463, groups: Clutch:Parental.tank, 16; Parental.tank, 14
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.540994 0.007851 14.012366 68.906 <2e-16 ***
## TreatmentEmbryo 0.002915 0.011338 11.344249 0.257 0.802
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## TrtmntEmbry -0.627anova
Anova(lmer.HatchSL, type = 3)## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: SL
## Chisq Df Pr(>Chisq)
## (Intercept) 4748.0720 1 <2e-16 ***
## Treatment 0.0661 1 0.7971
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Emmean
EmmHatchSL = (emmeans(lmer.HatchSL, ~ Treatment) %>% as.data.frame)
EmmHatchSL## Treatment emmean SE df lower.CL upper.CL
## Control 0.5409943 0.008693179 13.95 0.5223429 0.5596456
## Embryo 0.5439093 0.010178916 13.65 0.5220247 0.5657939
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.956. Figure
fig.HatchSL <- ggplot(EmmHatchSL,
aes(x = Treatment, y = emmean, fill = Treatment)) +
geom_errorbar(
aes(ymin = lower.CL, ymax = upper.CL),
position = position_dodge(width = 0.6),
width = 0.1, colour = "black", linewidth = 0.3
) +
geom_point(
shape = 22, size = 8, colour = "black", stroke = 0.3,
position = position_dodge(width = 0.6)
) +
theme_classic() +
labs(
x = NULL,
y = "Standard legnth (cm)"
) +
scale_fill_manual(values = c("Control" = "#0072B2",
"Embryo" = "#B2182B")) +
scale_y_continuous(
limits = c(0.5, 0.6),
breaks = c(0.5, 0.55, 0.60),
labels = c("0.50", "0.550", "0.60")) +
theme(
text = element_text(family = "Helvetica"),
axis.title = element_text(size = 14, face = "bold", colour = "black"),
axis.text = element_text(size = 12, colour = "black"),
axis.text.x = element_blank(),
legend.position = "none"
)
fig.HatchSL
Juvenile weight at hatching
1. Random effect structure
lmer.HatchW = lmer(W ~ Treatment + (1|Clutch), data = Hatching)
lmer.HatchWa = lmer(W ~ Treatment + (1|Parental.tank/Clutch), data = Hatching)## boundary (singular) fit: see help('isSingular')AIC(lmer.HatchW, lmer.HatchWa)## df AIC
## lmer.HatchW 4 -5981.906
## lmer.HatchWa 5 -5979.906As per above, choosing “a” to account for same clutches from the same pair.
2. Covariates
lm = lm(W ~ FemaleSL, data=Hatching)
summary(lm)##
## Call:
## lm(formula = W ~ FemaleSL, data = Hatching)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.151e-03 -4.484e-04 -7.963e-05 3.277e-04 2.360e-03
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.292e-03 2.967e-04 4.355 1.64e-05 ***
## FemaleSL 2.132e-04 3.133e-05 6.805 3.15e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0005994 on 461 degrees of freedom
## (6 observations deleted due to missingness)
## Multiple R-squared: 0.09128, Adjusted R-squared: 0.08931
## F-statistic: 46.31 on 1 and 461 DF, p-value: 3.151e-11lm = lm(W ~ FemaleFK, data=Hatching)
summary(lm)##
## Call:
## lm(formula = W ~ FemaleFK, data = Hatching)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.184e-03 -4.655e-04 -8.396e-05 3.421e-04 2.387e-03
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.474e-03 3.825e-04 6.467 2.55e-10 ***
## FemaleFK 1.661e-04 7.645e-05 2.172 0.0304 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0006256 on 461 degrees of freedom
## (6 observations deleted due to missingness)
## Multiple R-squared: 0.01013, Adjusted R-squared: 0.007982
## F-statistic: 4.718 on 1 and 461 DF, p-value: 0.03037lm = lm(W ~ MaleSL, data=Hatching)
summary(lm)##
## Call:
## lm(formula = W ~ MaleSL, data = Hatching)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0011565 -0.0004151 -0.0001031 0.0002936 0.0024232
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.067e-03 3.951e-04 2.700 0.00719 **
## MaleSL 2.358e-04 4.157e-05 5.671 2.51e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0006079 on 461 degrees of freedom
## (6 observations deleted due to missingness)
## Multiple R-squared: 0.06521, Adjusted R-squared: 0.06318
## F-statistic: 32.16 on 1 and 461 DF, p-value: 2.509e-08lm = lm(W ~ MaleFK, data=Hatching)
summary(lm)##
## Call:
## lm(formula = W ~ MaleFK, data = Hatching)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.280e-03 -4.334e-04 -9.417e-05 3.118e-04 2.406e-03
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.832e-03 3.853e-04 7.351 9.04e-13 ***
## MaleFK 9.320e-05 7.623e-05 1.223 0.222
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0006277 on 461 degrees of freedom
## (6 observations deleted due to missingness)
## Multiple R-squared: 0.003232, Adjusted R-squared: 0.00107
## F-statistic: 1.495 on 1 and 461 DF, p-value: 0.2221Strongest correlation with Female SL. Testing model fit.
lmer.HatchWa = lmer(W ~ Treatment + (1|Parental.tank/Clutch), data = Hatching)## boundary (singular) fit: see help('isSingular')lmer.HatchWb = lmer(W ~ Treatment + FemaleSL + (1|Parental.tank/Clutch), data = Hatching) # fails to converge## boundary (singular) fit: see help('isSingular')## Warning: Model failed to converge with 1 negative eigenvalue: -3.0e+01b fails to converge. Choosing a (without covariate).
3. Final model
lmer.HatchW = lmer(W ~ Treatment + (1|Parental.tank/Clutch), data = Hatching)## boundary (singular) fit: see help('isSingular')4. Model performance
performance::check_model(lmer.HatchW, check="homogeneity") 
performance::check_model(lmer.HatchW, check="outliers") 
performance::check_model(lmer.HatchW, check="qq", detrend = FALSE) 
performance::check_model(lmer.HatchW, check="normality") 
performance::check_model(lmer.HatchW, check="linearity") 
performance::check_model(lmer.HatchW, check="pp_check", re_formula=NA)
hist(residuals(lmer.HatchW), col="darkgray")
shapiro.test(residuals(lmer.HatchW))##
## Shapiro-Wilk normality test
##
## data: residuals(lmer.HatchW)
## W = 0.99749, p-value = 0.7145qqline(resid(lmer.HatchW))
outlierTest(lmer.HatchW)## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
## rstudent unadjusted p-value Bonferroni p
## 89 -3.255552 0.0012158 0.56292#89 identified as an outlier but retained as it appears within expected range of that clutch.
5. Output
Summary statistics
summary(lmer.HatchW)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: W ~ Treatment + (1 | Parental.tank/Clutch)
## Data: Hatching
##
## REML criterion at convergence: -5989.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1995 -0.6487 0.0082 0.6510 3.0215
##
## Random effects:
## Groups Name Variance Std.Dev.
## Clutch:Parental.tank (Intercept) 3.134e-07 0.0005598
## Parental.tank (Intercept) 0.000e+00 0.0000000
## Residual 1.139e-07 0.0003375
## Number of obs: 463, groups: Clutch:Parental.tank, 16; Parental.tank, 14
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.328e-03 1.878e-04 1.401e+01 17.718 5.44e-11 ***
## TreatmentEmbryo -8.469e-05 2.839e-04 1.400e+01 -0.298 0.77
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## TrtmntEmbry -0.662
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')anova
Anova(lmer.HatchW, type = 3)## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: W
## Chisq Df Pr(>Chisq)
## (Intercept) 313.940 1 <2e-16 ***
## Treatment 0.089 1 0.7655
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Emmean
EmmHatchW = (emmeans(lmer.HatchW, ~ Treatment) %>% as.data.frame)
EmmHatchW## Treatment emmean SE df lower.CL upper.CL
## Control 0.003327832 0.0002061064 14.01 0.002885802 0.003769861
## Embryo 0.003243141 0.0002413916 13.99 0.002725369 0.003760913
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.956. Figure
fig.HatchW <- ggplot(EmmHatchW, aes(x = Treatment, y = emmean, fill = Treatment)) +
geom_errorbar(
aes(ymin = lower.CL, ymax = upper.CL),
position = position_dodge(width = 0.6),
width = 0.1, colour = "black", linewidth = 0.3
) +
geom_point(
shape = 22, size = 8, colour = "black", stroke = 0.3,
position = position_dodge(width = 0.6)
) +
theme_classic() +
labs(
x = NULL,
y = "Weight (g)"
) +
scale_fill_manual(values = c("Control" = "#0072B2", "Embryo" = "#B2182B")) +
scale_y_continuous(
limits = c(0.002, 0.004),
breaks = c(0.002, 0.003, 0.004),
labels = c(0.002, 0.003, 0.004)) +
theme(
text = element_text(family = "Helvetica"),
axis.title.x = element_text(size = 14, face = "bold", colour = "black", margin = margin(t = 8)),
axis.title.y = element_text(size = 14, face = "bold", colour = "black", margin = margin(r = 8)),
axis.text = element_text(size = 12, colour = "black"),
axis.text.x = element_blank(),
legend.position = "none"
)
fig.HatchW
Condition (W for given SL) at hatching
1. Random effect structure
lmer.HatchCond = lmer(W ~ Treatment + SL + (1|Clutch), data = Hatching)
lmer.HatchConda = lmer(W ~ Treatment + SL + (1|Parental.tank/Clutch), data = Hatching)## boundary (singular) fit: see help('isSingular')AIC(lmer.HatchCond, lmer.HatchConda)## df AIC
## lmer.HatchCond 5 -5998.879
## lmer.HatchConda 6 -5996.879Choosing “a” to account for cases where clutches from the same pair is used.
2. Covariate exploration
Testing Female SL to follow hatching W.
lmer.HatchConda = lmer(W ~ Treatment + SL + (1|Parental.tank/Clutch), data = Hatching)## boundary (singular) fit: see help('isSingular')lmer.HatchCondb = lmer(W ~ Treatment + SL + FemaleSL + (1|Parental.tank/Clutch), data = Hatching)## boundary (singular) fit: see help('isSingular')AIC(lmer.HatchConda, lmer.HatchCondb)## df AIC
## lmer.HatchConda 6 -5996.879
## lmer.HatchCondb 7 -5981.132Choosing a (without covariate).
3. Final model
lmer.HatchCond = lmer(W ~ Treatment + SL + (1|Parental.tank/Clutch), data = Hatching)## boundary (singular) fit: see help('isSingular')4. Model performance
performance::check_model(lmer.HatchCond, check="homogeneity") 
performance::check_model(lmer.HatchCond, check="outliers") 
performance::check_model(lmer.HatchCond, check="qq", detrend = FALSE) 
performance::check_model(lmer.HatchCond, check="normality") 
performance::check_model(lmer.HatchCond, check="linearity") 
performance::check_model(lmer.HatchCond, check="pp_check", re_formula=NA)
hist(residuals(lmer.HatchCond), col="darkgray")
shapiro.test(residuals(lmer.HatchCond))##
## Shapiro-Wilk normality test
##
## data: residuals(lmer.HatchCond)
## W = 0.99823, p-value = 0.9207qqline(resid(lmer.HatchCond))
outlierTest(lmer.HatchCond)## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
## rstudent unadjusted p-value Bonferroni p
## 155 3.210153 0.0014203 0.65761#155 identified as an outlier but retained as it appears within the range for that clutch.
5. Output
Summary statistics
summary(lmer.HatchCond)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: W ~ Treatment + SL + (1 | Parental.tank/Clutch)
## Data: Hatching
##
## REML criterion at convergence: -6008.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.05529 -0.67038 0.01595 0.65933 3.15689
##
## Random effects:
## Groups Name Variance Std.Dev.
## Clutch:Parental.tank (Intercept) 2.311e-07 0.0004807
## Parental.tank (Intercept) 0.000e+00 0.0000000
## Residual 1.074e-07 0.0003277
## Number of obs: 463, groups: Clutch:Parental.tank, 16; Parental.tank, 14
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.048e-03 4.281e-04 2.942e+02 2.448 0.015 *
## TreatmentEmbryo -1.001e-04 2.442e-04 1.372e+01 -0.410 0.688
## SL 4.212e-03 7.324e-04 4.577e+02 5.751 1.63e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) TrtmnE
## TrtmntEmbry -0.240
## SL -0.926 -0.011
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')anova
Anova(lmer.HatchCond, type = 3)## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: W
## Chisq Df Pr(>Chisq)
## (Intercept) 5.9912 1 0.01438 *
## Treatment 0.1679 1 0.68200
## SL 33.0704 1 8.888e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Emmean
EmmHatchCond = (emmeans(lmer.HatchCond, ~ Treatment) %>% as.data.frame)
EmmHatchCond## Treatment emmean SE df lower.CL upper.CL
## Control 0.003336068 0.0001773706 13.98 0.002955585 0.003716551
## Embryo 0.003236013 0.0002076922 13.95 0.002790411 0.003681616
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.956. Figure
fig.HatchCond <- ggplot(EmmHatchCond,
aes(x = Treatment, y = emmean, fill = Treatment)) +
geom_errorbar(
aes(ymin = lower.CL, ymax = upper.CL),
position = position_dodge(width = 0.6),
width = 0.1, colour = "black", linewidth = 0.3
) +
geom_point(
shape = 22, size = 8, colour = "black", stroke = 0.3,
position = position_dodge(width = 0.6)
) +
theme_classic() +
labs(
y = "Condition (g)"
) +
scale_y_continuous(
limits = c(0.002, 0.004),
breaks = c(0.002, 0.003, 0.004),
labels = c(0.002, 0.003, 0.004)) +
scale_fill_manual(values = c("Control" = "#0072B2",
"Embryo" = "#B2182B")) +
theme(
text = element_text(family = "Helvetica"),
axis.title.x = element_blank(),
axis.title = element_text(size = 14, face = "bold", colour = "black"),
axis.text = element_text(size = 12, colour = "black"),
legend.position = "none"
)
fig.HatchCond
Yolk area at hatching
1. Random effect structure
lmer.HatchYA = lmer(YA ~ Treatment + (1|Clutch), data = Hatching)
lmer.HatchYAa = lmer(YA ~ Treatment + (1|Parental.tank/Clutch), data = Hatching)
AIC(lmer.HatchYA, lmer.HatchYAa)## df AIC
## lmer.HatchYA 4 -4327.780
## lmer.HatchYAa 5 -4327.068Choosing “a” to account for cases where clutches from the same pair is used.
2. Covariates
lm = lm(YA ~ FemaleSL, data=Hatching)
summary(lm)##
## Call:
## lm(formula = YA ~ FemaleSL, data = Hatching)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.008650 -0.002663 -0.000203 0.002639 0.009384
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.0010061 0.0017774 -0.566 0.572
## FemaleSL 0.0016669 0.0001876 8.886 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.003585 on 465 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.1452, Adjusted R-squared: 0.1433
## F-statistic: 78.96 on 1 and 465 DF, p-value: < 2.2e-16lm = lm(YA ~ FemaleFK, data=Hatching)
summary(lm)##
## Call:
## lm(formula = YA ~ FemaleFK, data = Hatching)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0094994 -0.0027540 0.0004418 0.0026997 0.0096518
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0064127 0.0023222 2.761 0.005982 **
## FemaleFK 0.0016637 0.0004637 3.588 0.000369 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.003825 on 465 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.02693, Adjusted R-squared: 0.02484
## F-statistic: 12.87 on 1 and 465 DF, p-value: 0.0003692lm = lm(YA ~ MaleSL, data=Hatching)
summary(lm)##
## Call:
## lm(formula = YA ~ MaleSL, data = Hatching)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0093483 -0.0028937 -0.0001719 0.0022978 0.0102887
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0030975 0.0024560 1.261 0.208
## MaleSL 0.0012251 0.0002582 4.744 2.79e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.003787 on 465 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.04617, Adjusted R-squared: 0.04412
## F-statistic: 22.51 on 1 and 465 DF, p-value: 2.789e-06lm = lm(YA ~ MaleFK, data=Hatching)
summary(lm)##
## Call:
## lm(formula = YA ~ MaleFK, data = Hatching)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0086351 -0.0032779 -0.0002471 0.0027458 0.0096113
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.004761 0.002310 2.061 0.0398 *
## MaleFK 0.001976 0.000457 4.324 1.88e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.003802 on 465 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.03865, Adjusted R-squared: 0.03659
## F-statistic: 18.7 on 1 and 465 DF, p-value: 1.876e-05Strongest correlation with Female SL. Testing model fit.
lmer.HatchYAa = lmer(YA ~ Treatment + (1|Parental.tank/Clutch), data = Hatching)
lmer.HatchYAb = lmer(YA ~ Treatment + FemaleSL + (1|Parental.tank/Clutch), data = Hatching)
AIC(lmer.HatchYAa, lmer.HatchYAb)## df AIC
## lmer.HatchYAa 5 -4327.068
## lmer.HatchYAb 6 -4316.111Choosing a (without covariate).
3. Final model
lmer.HatchYA = lmer(YA ~ Treatment + (1|Parental.tank/Clutch), data = Hatching)4. Model performance
performance::check_model(lmer.HatchYA, check="homogeneity") 
performance::check_model(lmer.HatchYA, check="outliers") 
performance::check_model(lmer.HatchYA, check="qq", detrend = FALSE) 
performance::check_model(lmer.HatchYA, check="normality") 
performance::check_model(lmer.HatchYA, check="linearity") 
performance::check_model(lmer.HatchYA, check="pp_check", re_formula=NA)
hist(residuals(lmer.HatchYA), col="darkgray")
shapiro.test(residuals(lmer.HatchYA))##
## Shapiro-Wilk normality test
##
## data: residuals(lmer.HatchYA)
## W = 0.99535, p-value = 0.1779qqline(resid(lmer.HatchYA))
outlierTest(lmer.HatchYA)## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
## rstudent unadjusted p-value Bonferroni p
## 345 3.165013 0.0016534 0.77215#345 identified as an outlier but retained as it appears within the expected range of the clutch.
5. Model outputs
Summary statistics
summary(lmer.HatchYA)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: YA ~ Treatment + (1 | Parental.tank/Clutch)
## Data: Hatching
##
## REML criterion at convergence: -4337.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.03325 -0.57097 0.07431 0.61691 3.11098
##
## Random effects:
## Groups Name Variance Std.Dev.
## Clutch:Parental.tank (Intercept) 3.909e-06 0.001977
## Parental.tank (Intercept) 7.503e-06 0.002739
## Residual 4.477e-06 0.002116
## Number of obs: 467, groups: Clutch:Parental.tank, 16; Parental.tank, 14
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.0145352 0.0010832 13.6866462 13.419 2.91e-09 ***
## TreatmentEmbryo -0.0004212 0.0014146 5.9970262 -0.298 0.776
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## TrtmntEmbry -0.563anova
Anova(lmer.HatchYA, type = 3)## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: YA
## Chisq Df Pr(>Chisq)
## (Intercept) 180.0662 1 <2e-16 ***
## Treatment 0.0887 1 0.7659
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Emmean
EmmHatchYA = (emmeans(lmer.HatchYA, ~ Treatment) %>% as.data.frame)
EmmHatchYA## Treatment emmean SE df lower.CL upper.CL
## Control 0.01453521 0.001216352 13.54 0.01191809 0.01715233
## Embryo 0.01411403 0.001414322 12.13 0.01103622 0.01719184
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.956. Figure
fig.HatchYA <- ggplot(EmmHatchYA,
aes(x = Treatment, y = emmean, fill = Treatment)) +
geom_errorbar(
aes(ymin = lower.CL, ymax = upper.CL),
position = position_dodge(width = 0.6),
width = 0.1, colour = "black", linewidth = 0.3
) +
geom_point(
shape = 22, size = 8, colour = "black", stroke = 0.3,
position = position_dodge(width = 0.6)
) +
theme_classic() +
labs(
x = NULL,
y = "Yolk area (cm²)"
) +
scale_fill_manual(values = c("Control" = "#0072B2",
"Embryo" = "#B2182B")) +
scale_y_continuous(
limits = c(0.01, 0.02),
breaks = c(0.01, 0.015,0.02),
labels = c(0.01, 0.015,0.02)) +
theme(
text = element_text(family = "Helvetica"),
axis.title = element_text(size = 14, face = "bold", colour = "black"),
axis.text = element_text(size = 12, colour = "black"),
legend.position = "none"
)
fig.HatchYA
Below are effects of MHWs during embryonic development on post-hatching growth.
Growth in standard length
1. Random effect structures
lmer.SLa = lmer(LogSL ~ Treatment * Day + (1|Tank) + (1|Parental.tank/Clutch), data = MHWembryo)## boundary (singular) fit: see help('isSingular')lmer.SLb = lmer(LogSL ~ Treatment * Day + (1|Parental.tank/Clutch/Tank), data = MHWembryo)## boundary (singular) fit: see help('isSingular')AIC(lmer.SLa, lmer.SLb)## df AIC
## lmer.SLa 12 -2256.799
## lmer.SLb 12 -2275.611lmer.SL and lmer.SLb essentially a same model. Opting for the nested model.
2. Covariate exploration
lm = lm(SL ~ FemaleSL, data=MHWembryo)
summary(lm)##
## Call:
## lm(formula = SL ~ FemaleSL, data = MHWembryo)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.90743 -0.54934 -0.03995 0.46605 1.17393
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.485570 0.160678 9.246 <2e-16 ***
## FemaleSL 0.000725 0.017098 0.042 0.966
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5312 on 1159 degrees of freedom
## (21 observations deleted due to missingness)
## Multiple R-squared: 1.551e-06, Adjusted R-squared: -0.0008613
## F-statistic: 0.001798 on 1 and 1159 DF, p-value: 0.9662lm = lm(SL ~ FemaleFK, data=MHWembryo)
summary(lm)##
## Call:
## lm(formula = SL ~ FemaleFK, data = MHWembryo)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.90954 -0.54018 -0.04186 0.46475 1.20085
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.81816 0.20443 8.894 <2e-16 ***
## FemaleFK -0.06589 0.04122 -1.598 0.11
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5306 on 1159 degrees of freedom
## (21 observations deleted due to missingness)
## Multiple R-squared: 0.0022, Adjusted R-squared: 0.001339
## F-statistic: 2.555 on 1 and 1159 DF, p-value: 0.1102lm = lm(SL ~ MaleSL, data=MHWembryo)
summary(lm)##
## Call:
## lm(formula = SL ~ MaleSL, data = MHWembryo)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.90682 -0.54926 -0.03959 0.46641 1.17417
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.4876067 0.2134145 6.971 5.29e-12 ***
## MaleSL 0.0005031 0.0225733 0.022 0.982
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5312 on 1159 degrees of freedom
## (21 observations deleted due to missingness)
## Multiple R-squared: 4.286e-07, Adjusted R-squared: -0.0008624
## F-statistic: 0.0004968 on 1 and 1159 DF, p-value: 0.9822lm = lm(SL ~ MaleFK, data=MHWembryo)
summary(lm)##
## Call:
## lm(formula = SL ~ MaleFK, data = MHWembryo)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.88774 -0.54287 -0.02382 0.47302 1.18526
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.68163 0.22009 7.641 4.5e-14 ***
## MaleFK -0.03763 0.04365 -0.862 0.389
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.531 on 1159 degrees of freedom
## (21 observations deleted due to missingness)
## Multiple R-squared: 0.000641, Adjusted R-squared: -0.0002212
## F-statistic: 0.7434 on 1 and 1159 DF, p-value: 0.3887No parental trait with significant correlation.
3. Final model
lmer.SL = lmer(LogSL ~ Treatment * Day + (1|Parental.tank/Clutch/Tank), data = MHWembryo)## boundary (singular) fit: see help('isSingular')4. 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", re_formula=NA)
hist(residuals(lmer.SL), col="darkgray")
shapiro.test(residuals(lmer.SL))##
## Shapiro-Wilk normality test
##
## data: residuals(lmer.SL)
## W = 0.98902, p-value = 1.235e-07qqline(resid(lmer.SL))
outlierTest(lmer.SL)## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
## rstudent unadjusted p-value Bonferroni p
## 552 -4.010776 6.4432e-05 0.074805#522 identified as outlier but retained as it appears within the expected range for that clutch.
5. Model outouts
Summary statistics
summary(lmer.SL)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: LogSL ~ Treatment * Day + (1 | Parental.tank/Clutch/Tank)
## Data: MHWembryo
##
## REML criterion at convergence: -2299.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.9048 -0.5574 0.0191 0.6235 3.4499
##
## Random effects:
## Groups Name Variance Std.Dev.
## Tank:(Clutch:Parental.tank) (Intercept) 4.980e-03 7.057e-02
## Clutch:Parental.tank (Intercept) 2.981e-11 5.460e-06
## Parental.tank (Intercept) 8.030e-11 8.961e-06
## Residual 6.766e-03 8.225e-02
## Number of obs: 1161, groups:
## Tank:(Clutch:Parental.tank), 64; Clutch:Parental.tank, 17; Parental.tank, 15
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -0.160149 0.024455 56.698612 -6.549 1.83e-08 ***
## TreatmentEmbryo -0.079528 0.035597 56.379756 -2.234 0.0295 *
## Day30 0.428217 0.034520 56.283745 12.405 < 2e-16 ***
## Day45 0.773105 0.035635 56.617655 21.695 < 2e-16 ***
## Day60 0.926486 0.035668 56.822672 25.975 < 2e-16 ***
## TreatmentEmbryo:Day30 0.050265 0.050276 56.088879 1.000 0.3217
## TreatmentEmbryo:Day45 0.003362 0.053201 56.297294 0.063 0.9498
## TreatmentEmbryo:Day60 0.063104 0.051068 56.339418 1.236 0.2217
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) TrtmnE Day30 Day45 Day60 TE:D30 TE:D45
## TrtmntEmbry -0.687
## Day30 -0.708 0.487
## Day45 -0.686 0.471 0.486
## Day60 -0.686 0.471 0.486 0.471
## TrtmntE:D30 0.486 -0.708 -0.687 -0.334 -0.334
## TrtmntE:D45 0.460 -0.669 -0.326 -0.670 -0.315 0.474
## TrtmntE:D60 0.479 -0.697 -0.339 -0.329 -0.698 0.494 0.466
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')anova
Anova(lmer.SL, type = 3)## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: LogSL
## Chisq Df Pr(>Chisq)
## (Intercept) 42.8850 1 5.805e-11 ***
## Treatment 4.9912 1 0.02548 *
## Day 806.4622 3 < 2.2e-16 ***
## Treatment:Day 2.3074 3 0.51111
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Emmean
EmmSL.Log = (emmeans(lmer.SL, ~ Treatment * Day) %>% as.data.frame)
# Back-transforming
EmmSL <- EmmSL.Log %>%
mutate(
emmean = exp(emmean),
lower.CL = exp(lower.CL),
upper.CL = exp(upper.CL)
)
EmmSL## Treatment Day emmean SE df lower.CL upper.CL
## 1 Control 15 0.8520171 0.02517301 56.49105 0.8101251 0.8960754
## 2 Embryo 15 0.7868819 0.02667223 55.89152 0.7459399 0.8300711
## 3 Control 30 1.3074360 0.02508896 55.66558 1.2433409 1.3748351
## 4 Embryo 30 1.2697309 0.02662463 55.53091 1.2037715 1.3393045
## 5 Control 45 1.8458799 0.02653703 56.33871 1.7503278 1.9466483
## 6 Embryo 45 1.7105074 0.03136071 55.78965 1.6063455 1.8214235
## 7 Control 60 2.1518707 0.02690709 56.72470 2.0389835 2.2710079
## 8 Embryo 60 2.1168163 0.02661693 55.46887 2.0068811 2.2327736Pairwise comparison
Day
emm_day <- emmeans(lmer.SL, ~ Day)## NOTE: Results may be misleading due to involvement in interactionspairs_day <- pairs(emm_day, adjust = "tukey") |> summary(infer = TRUE)
pairs_day## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Day15 - Day30 -0.453 0.0251 41.1 -0.521 -0.386 -18.034 <.0001
## Day15 - Day45 -0.775 0.0268 43.6 -0.846 -0.703 -28.870 <.0001
## Day15 - Day60 -0.958 0.0256 41.9 -1.026 -0.890 -37.475 <.0001
## Day30 - Day45 -0.321 0.0268 43.4 -0.393 -0.250 -11.992 <.0001
## Day30 - Day60 -0.505 0.0255 41.7 -0.573 -0.436 -19.767 <.0001
## Day45 - Day60 -0.183 0.0272 44.2 -0.256 -0.111 -6.730 <.0001
##
## Results are averaged over the levels of: Treatment
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Conf-level adjustment: tukey method for comparing a family of 4 estimates
## P value adjustment: tukey method for comparing a family of 4 estimates2. Figure
EmmSL$Day <- as.numeric(as.character(EmmSL$Day))
fig.Embryo.SL <- ggplot(EmmSL,
aes(x = Day, y = emmean, colour = Treatment, fill = Treatment)) +
geom_ribbon(aes(ymin = lower.CL, ymax = upper.CL),
alpha = 0.15, colour = NA) +
geom_line(linewidth = 1) +
geom_point(size = 3, stroke = 0.5, colour = "black",
aes(fill = Treatment), shape = 22) +
scale_fill_manual(values = c("Control" = "#0072B2",
"Embryo" = "#B2182B")) +
scale_colour_manual(values = c("Control" = "#0072B2",
"Embryo" = "#B2182B")) +
scale_x_continuous(breaks = c(15, 30, 45, 60),
labels = c(15, 30, 45, 60)) +
scale_y_continuous(
limits = c(0.5, 2.5),
breaks = c(0.5, 1, 1.5, 2, 2.5),
labels = c(0.5, 1, 1.5, 2, 2.5)) +
theme_classic() +
labs(
y = "Standard length (cm)",
x = NULL
) +
theme(
text = element_text(family = "Helvetica"),
axis.title = element_text(size = 14, face = "bold", colour = "black"),
axis.text = element_text(size = 12, colour = "black"),
axis.text.x = element_blank(),
legend.position = "none"
)
fig.Embryo.SL
Growth in Weight
1. Random effect structure
lmer.Wa = lmer(LogW ~ Treatment * Day + (1|Tank) + (1|Parental.tank/Clutch), data = MHWembryo)## boundary (singular) fit: see help('isSingular')lmer.Wb = lmer(LogW ~ Treatment * Day + (1|Parental.tank/Clutch/Tank), data = MHWembryo)
AIC(lmer.Wa, lmer.Wb)## df AIC
## lmer.Wa 12 412.9182
## lmer.Wb 12 363.9054Choosing b.
2. Covariate exploration
lm = lm(W ~ FemaleSL, data=MHWembryo)
summary(lm)##
## Call:
## lm(formula = W ~ FemaleSL, data = MHWembryo)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.16535 -0.12963 -0.05558 0.10911 0.51875
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.115643 0.045315 2.552 0.0108 *
## FemaleSL 0.005322 0.004822 1.104 0.2699
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1497 on 1156 degrees of freedom
## (24 observations deleted due to missingness)
## Multiple R-squared: 0.001053, Adjusted R-squared: 0.0001888
## F-statistic: 1.219 on 1 and 1156 DF, p-value: 0.2699lm = lm(W ~ FemaleFK, data=MHWembryo)
summary(lm)##
## Call:
## lm(formula = W ~ FemaleFK, data = MHWembryo)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.16012 -0.12992 -0.05623 0.10709 0.52610
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.22388 0.05778 3.875 0.000113 ***
## FemaleFK -0.01182 0.01165 -1.015 0.310507
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1497 on 1156 degrees of freedom
## (24 observations deleted due to missingness)
## Multiple R-squared: 0.0008897, Adjusted R-squared: 2.543e-05
## F-statistic: 1.029 on 1 and 1156 DF, p-value: 0.3105lm = lm(W ~ MaleSL, data=MHWembryo)
summary(lm)##
## Call:
## lm(formula = W ~ MaleSL, data = MHWembryo)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.16578 -0.12896 -0.05423 0.10804 0.51832
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.081094 0.060145 1.348 0.178
## MaleSL 0.008944 0.006362 1.406 0.160
##
## Residual standard error: 0.1496 on 1156 degrees of freedom
## (24 observations deleted due to missingness)
## Multiple R-squared: 0.001707, Adjusted R-squared: 0.0008435
## F-statistic: 1.977 on 1 and 1156 DF, p-value: 0.16lm = lm(W ~ MaleFK, data=MHWembryo)
summary(lm)##
## Call:
## lm(formula = W ~ MaleFK, data = MHWembryo)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.15669 -0.12998 -0.05411 0.10991 0.53171
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.21675 0.06206 3.493 0.000496 ***
## MaleFK -0.01020 0.01231 -0.829 0.407203
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1497 on 1156 degrees of freedom
## (24 observations deleted due to missingness)
## Multiple R-squared: 0.0005943, Adjusted R-squared: -0.0002702
## F-statistic: 0.6875 on 1 and 1156 DF, p-value: 0.4072No parental trait with significant correlation.
3. Final model
lmer.W = lmer(LogW ~ Treatment * Day + (1|Parental.tank/Clutch/Tank), data = MHWembryo)4. 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")## Failed to compute posterior predictive checks with `re_formula=NULL`.
## Trying again with `re_formula=NA` now.
hist(residuals(lmer.W), col="darkgray")
shapiro.test(residuals(lmer.W))##
## Shapiro-Wilk normality test
##
## data: residuals(lmer.W)
## W = 0.98988, p-value = 3.708e-07qqline(resid(lmer.W))
outlierTest(lmer.W)## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
## rstudent unadjusted p-value Bonferroni p
## 552 -3.958992 7.9906e-05 0.092531#552 identified as an outlier but retained as it appears within the range for that clutch.
5. Model output
Summary statistics
summary(lmer.W)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: LogW ~ Treatment * Day + (1 | Parental.tank/Clutch/Tank)
## Data: MHWembryo
##
## REML criterion at convergence: 339.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.8528 -0.5593 -0.0081 0.6352 3.1428
##
## Random effects:
## Groups Name Variance Std.Dev.
## Tank:(Clutch:Parental.tank) (Intercept) 6.374e-02 2.525e-01
## Clutch:Parental.tank (Intercept) 4.029e-09 6.347e-05
## Parental.tank (Intercept) 2.199e-09 4.689e-05
## Residual 6.597e-02 2.568e-01
## Number of obs: 1158, groups:
## Tank:(Clutch:Parental.tank), 64; Clutch:Parental.tank, 17; Parental.tank, 15
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -3.82366 0.08674 56.63947 -44.082 < 2e-16 ***
## TreatmentEmbryo -0.09690 0.12630 56.38105 -0.767 0.446
## Day30 1.24532 0.12247 56.27609 10.169 2.35e-14 ***
## Day45 2.37771 0.12638 56.53229 18.814 < 2e-16 ***
## Day60 2.78203 0.12647 56.69365 21.997 < 2e-16 ***
## TreatmentEmbryo:Day30 0.08357 0.17840 56.11943 0.468 0.641
## TreatmentEmbryo:Day45 -0.21529 0.18873 56.27720 -1.141 0.259
## TreatmentEmbryo:Day60 0.04231 0.18116 56.31370 0.234 0.816
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) TrtmnE Day30 Day45 Day60 TE:D30 TE:D45
## TrtmntEmbry -0.687
## Day30 -0.708 0.486
## Day45 -0.686 0.471 0.486
## Day60 -0.686 0.471 0.486 0.471
## TrtmntE:D30 0.486 -0.708 -0.686 -0.334 -0.333
## TrtmntE:D45 0.460 -0.669 -0.326 -0.670 -0.315 0.474
## TrtmntE:D60 0.479 -0.697 -0.339 -0.329 -0.698 0.494 0.467anova
Anova(lmer.W, type = 3)## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: LogW
## Chisq Df Pr(>Chisq)
## (Intercept) 1943.2143 1 <2e-16 ***
## Treatment 0.5887 1 0.4429
## Day 591.9820 3 <2e-16 ***
## Treatment:Day 2.8598 3 0.4137
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Pairwise comparison
Day
emm_day <- emmeans(lmer.W, ~ Day)## NOTE: Results may be misleading due to involvement in interactionspairs_day <- pairs(emm_day, adjust = "tukey") |> summary(infer = TRUE)
pairs_day## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Day15 - Day30 -1.287 0.0892 41.1 -1.526 -1.048 -14.429 <.0001
## Day15 - Day45 -2.270 0.0952 43.5 -2.524 -2.016 -23.844 <.0001
## Day15 - Day60 -2.803 0.0907 41.9 -3.046 -2.561 -30.909 <.0001
## Day30 - Day45 -0.983 0.0951 43.4 -1.237 -0.729 -10.336 <.0001
## Day30 - Day60 -1.516 0.0906 41.7 -1.758 -1.274 -16.736 <.0001
## Day45 - Day60 -0.533 0.0966 44.1 -0.791 -0.275 -5.520 <.0001
##
## Results are averaged over the levels of: Treatment
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Conf-level adjustment: tukey method for comparing a family of 4 estimates
## P value adjustment: tukey method for comparing a family of 4 estimatesEmmean
EmmW.Log <- emmeans(lmer.W, ~ Treatment * Day) %>%
as.data.frame()
EmmW <- EmmW.Log %>%
mutate(
emmean = exp(emmean),
lower.CL = exp(lower.CL),
upper.CL = exp(upper.CL)
)
EmmW## Treatment Day emmean SE df lower.CL upper.CL
## 1 Control 15 0.02184776 0.08928752 56.44429 0.01827008 0.02612602
## 2 Embryo 15 0.01983004 0.09465566 55.95822 0.01640484 0.02397038
## 3 Control 30 0.07590006 0.08902353 55.72101 0.06350142 0.09071954
## 4 Embryo 30 0.07489522 0.09449623 55.61669 0.06197695 0.09050613
## 5 Control 45 0.23552255 0.09410678 56.24259 0.19505978 0.28437883
## 6 Embryo 45 0.17236492 0.11127268 55.81748 0.13792248 0.21540843
## 7 Control 60 0.35287888 0.09539812 56.54634 0.29150589 0.42717319
## 8 Embryo 60 0.33413301 0.09447493 55.56812 0.27651102 0.403762836. Figure
EmmW$Day <- as.numeric(as.character(EmmW$Day))
fig.Embryo.W <- ggplot(EmmW,
aes(x = Day, y = emmean, colour = Treatment, fill = Treatment)) +
geom_ribbon(aes(ymin = lower.CL, ymax = upper.CL),
alpha = 0.15, colour = NA) +
geom_line(linewidth = 1) +
geom_point(size = 3, stroke = 0.5, colour = "black",
aes(fill = Treatment), shape = 22) +
scale_fill_manual(values = c("Control" = "#0072B2",
"Embryo" = "#B2182B")) +
scale_colour_manual(values = c("Control" = "#0072B2",
"Embryo" = "#B2182B")) +
scale_x_continuous(breaks = c(15, 30, 45, 60),
labels = c(15, 30, 45, 60)) +
scale_y_continuous(
limits = c(0, 0.45),
breaks = c(0, 0.1, 0.2, 0.3, 0.4),
labels = c(0, 0.1, 0.2, 0.3, 0.4)) +
theme_classic() +
labs(
y = "Weight (g)",
x = NULL
) +
theme(
text = element_text(family = "Helvetica"),
axis.title = element_text(size = 14, face = "bold", colour = "black"),
axis.text = element_text(size = 12, colour = "black"),
axis.text.x = element_blank(),
legend.position = "none"
)
fig.Embryo.W
Growth in Condition (W for given SL)
MHWembryo <- MHWembryo %>%
group_by(Day) %>%
mutate(LogSL.centered = LogSL - mean(LogSL, na.rm = TRUE))1. Random effect structure
lmer.Conda = lmer(LogW ~ Treatment * Day + LogSL.centered + (1|Tank) + (1|Parental.tank/Clutch), data = MHWembryo) # fails to converge## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00331334 (tol = 0.002, component 1)lmer.Condb = lmer(LogW ~ Treatment * Day + LogSL.centered + (1|Parental.tank/Clutch/Tank), data = MHWembryo) # fails to converge## boundary (singular) fit: see help('isSingular')## Warning: Model failed to converge with 1 negative eigenvalue: -4.2e+01lmer.Condc = lmer(LogW ~ Treatment * Day + LogSL.centered + (1|Clutch/Tank), data = MHWembryo) ## boundary (singular) fit: see help('isSingular')a and b fails to converge. Choosing c.
2. Covariate explotation
No parental trait with significant correlation for W.
3. Final model
lmer.Cond = lmer(LogW ~ Treatment * Day + LogSL.centered + (1|Clutch/Tank), data = MHWembryo) ## boundary (singular) fit: see help('isSingular')4. Model performance
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")## Failed to compute posterior predictive checks with `re_formula=NULL`.
## Trying again with `re_formula=NA` now.
hist(residuals(lmer.Cond), col="darkgray")
shapiro.test(residuals(lmer.Cond))##
## Shapiro-Wilk normality test
##
## data: residuals(lmer.Cond)
## W = 0.97123, p-value = 2.117e-14qqline(resid(lmer.Cond))
outlierTest(lmer.Cond)## rstudent unadjusted p-value Bonferroni p
## 839 5.167133 2.8025e-07 0.00032453
## 1091 4.795175 1.8393e-06 0.00213000
## 473 4.621031 4.2514e-06 0.00492310
## 805 -4.506166 7.2777e-06 0.00842760
## 812 4.178724 3.1545e-05 0.03652900
## 843 4.139144 3.7416e-05 0.04332700All outliers appeared to be within the expected range for a particular clutch. Kept.
5. Model outputs
Summary statistics
summary(lmer.Cond)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: LogW ~ Treatment * Day + LogSL.centered + (1 | Clutch/Tank)
## Data: MHWembryo
##
## REML criterion at convergence: -2154.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.3886 -0.5961 -0.0004 0.5503 4.9848
##
## Random effects:
## Groups Name Variance Std.Dev.
## Tank:Clutch (Intercept) 0.014888 0.12202
## Clutch (Intercept) 0.000000 0.00000
## Residual 0.007255 0.08517
## Number of obs: 1158, groups: Tank:Clutch, 64; Clutch, 17
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -3.94148 0.04128 56.25568 -95.471 <2e-16 ***
## TreatmentEmbryo 0.13693 0.06017 56.21316 2.276 0.0267 *
## Day30 1.33165 0.05832 56.00480 22.834 <2e-16 ***
## Day45 2.39727 0.06014 56.10311 39.861 <2e-16 ***
## Day60 2.87072 0.06017 56.20838 47.711 <2e-16 ***
## LogSL.centered 2.95124 0.03103 1127.24521 95.113 <2e-16 ***
## TreatmentEmbryo:Day30 -0.06384 0.08499 55.93794 -0.751 0.4557
## TreatmentEmbryo:Day45 -0.22418 0.08986 55.97767 -2.495 0.0156 *
## TreatmentEmbryo:Day60 -0.14294 0.08627 56.05224 -1.657 0.1032
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) TrtmnE Day30 Day45 Day60 LgSL.c TE:D30 TE:D45
## TrtmntEmbry -0.687
## Day30 -0.708 0.486
## Day45 -0.686 0.471 0.486
## Day60 -0.686 0.471 0.486 0.471
## LogSL.cntrd -0.030 0.041 0.016 0.003 0.015
## TrtmntE:D30 0.486 -0.708 -0.686 -0.333 -0.333 -0.018
## TrtmntE:D45 0.459 -0.668 -0.325 -0.669 -0.315 -0.001 0.473
## TrtmntE:D60 0.479 -0.697 -0.339 -0.328 -0.698 -0.023 0.493 0.466
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')anova
Anova(lmer.Cond, type = 3)## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: LogW
## Chisq Df Pr(>Chisq)
## (Intercept) 9114.7372 1 < 2e-16 ***
## Treatment 5.1787 1 0.02287 *
## Day 2719.9320 3 < 2e-16 ***
## LogSL.centered 9046.5271 1 < 2e-16 ***
## Treatment:Day 7.0659 3 0.06983 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Emmean
EmmCond.Log <- emmeans(lmer.Cond, ~ Treatment * Day) |> as.data.frame()
EmmCond <- EmmCond.Log |>
mutate(
emmean = exp(emmean),
lower.CL = exp(lower.CL),
upper.CL = exp(upper.CL)
)
EmmCond## Treatment Day emmean SE df lower.CL upper.CL
## 1 Control 15 0.01943232 0.04128445 56.28062 0.01789004 0.02110756
## 2 Embryo 15 0.02228396 0.04373907 56.02881 0.02041456 0.02432455
## 3 Control 30 0.07359566 0.04120109 55.83311 0.06776494 0.07992808
## 4 Embryo 30 0.07917584 0.04369304 55.79510 0.07253988 0.08641886
## 5 Control 45 0.21362255 0.04394695 56.13975 0.19562101 0.23328063
## 6 Embryo 45 0.19577384 0.05097068 55.95046 0.17677038 0.21682025
## 7 Control 60 0.34297570 0.04396499 56.23570 0.31406353 0.37454948
## 8 Embryo 60 0.34092186 0.04368490 55.75428 0.31235284 0.37210393Pariwise
Day
emm_Cond <- emmeans(lmer.Cond, ~ Day)## NOTE: Results may be misleading due to involvement in interactionspairs_Cond <- pairs(emm_Cond, adjust = "tukey") |> summary(infer = TRUE)
pairs_Cond## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Day15 - Day30 -1.300 0.0425 41.1 -1.413 -1.186 -30.591 <.0001
## Day15 - Day45 -2.285 0.0451 43.5 -2.406 -2.165 -50.654 <.0001
## Day15 - Day60 -2.799 0.0432 41.9 -2.915 -2.684 -64.838 <.0001
## Day30 - Day45 -0.985 0.0451 43.4 -1.106 -0.865 -21.855 <.0001
## Day30 - Day60 -1.500 0.0431 41.8 -1.615 -1.384 -34.752 <.0001
## Day45 - Day60 -0.514 0.0457 44.1 -0.636 -0.392 -11.238 <.0001
##
## Results are averaged over the levels of: Treatment
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Conf-level adjustment: tukey method for comparing a family of 4 estimates
## P value adjustment: tukey method for comparing a family of 4 estimatesTreatment x Day
emm.Cond <- emmeans(lmer.Cond, ~ Treatment | Day)
pairs.Cond <- pairs(emm.Cond, adjust = "tukey") %>%
summary(infer = TRUE)
pairs.Cond## Day = 15:
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Control - Embryo -0.13693 0.0602 56.2 -0.2575 -0.0164 -2.276 0.0267
##
## Day = 30:
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Control - Embryo -0.07309 0.0601 55.8 -0.1934 0.0472 -1.217 0.2288
##
## Day = 45:
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Control - Embryo 0.08725 0.0673 56.1 -0.0476 0.2221 1.296 0.2003
##
## Day = 60:
## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Control - Embryo 0.00601 0.0620 56.0 -0.1182 0.1302 0.097 0.9231
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.956. Figure
EmmCond$Day <- as.numeric(as.character(EmmCond$Day))
fig.Embryo.Cond <- ggplot(EmmCond,
aes(x = Day, y = emmean, colour = Treatment, fill = Treatment)) +
geom_ribbon(aes(ymin = lower.CL, ymax = upper.CL),
alpha = 0.15, colour = NA) +
geom_line(linewidth = 1) +
geom_point(size = 3, stroke = 0.5, colour = "black",
aes(fill = Treatment), shape = 22) +
scale_fill_manual(values = c("Control" = "#0072B2",
"Embryo" = "#B2182B")) +
scale_colour_manual(values = c("Control" = "#0072B2",
"Embryo" = "#B2182B")) +
scale_x_continuous(breaks = c(15, 30, 45, 60),
labels = c(15, 30, 45, 60)) +
scale_y_continuous(
limits = c(0, 0.4),
breaks = c(0, 0.1, 0.2, 0.3, 0.4),
labels = c(0, 0.1, 0.2, 0.3, 0.4)) +
theme_classic() +
labs(
y = "Condition (g)",
x = NULL
) +
theme(
text = element_text(family = "Helvetica"),
axis.title = element_text(size = 14, face = "bold", colour = "black"),
axis.text = element_text(size = 12, colour = "black"),
legend.position = "none"
)
fig.Embryo.Cond
Effects of post-hatching MHWs on juvenile growth
Grwoth in SL (Day 15 and 30)
1. Random effect structures
Each clutch came from unique parents unlike impacts of embryonic MHWs above. Parental.tank are not included in the models below for this reason.
lmer.SL15a = lmer(LogSL ~ Treatment * Day + (1|Tank) + (1|Clutch), data = MHWday15_30)
lmer.SL15b = lmer(LogSL ~ Treatment * Day + (1|Clutch/Tank), data = MHWday15_30)
AIC(lmer.SL15a, lmer.SL15b)## df AIC
## lmer.SL15a 7 -954.5558
## lmer.SL15b 7 -958.6867Choosing b.
2. Covariate exploration
lm = lm(SL ~ FemaleSL, data=MHWday15_30)
summary(lm)##
## Call:
## lm(formula = SL ~ FemaleSL, data = MHWday15_30)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.36181 -0.20687 -0.08072 0.19150 0.56038
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.21519 0.42856 -2.836 0.00479 **
## FemaleSL 0.23796 0.04477 5.315 1.71e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2372 on 437 degrees of freedom
## (5 observations deleted due to missingness)
## Multiple R-squared: 0.06071, Adjusted R-squared: 0.05856
## F-statistic: 28.25 on 1 and 437 DF, p-value: 1.706e-07lm = lm(SL ~ FemaleFK, data=MHWday15_30)
summary(lm)##
## Call:
## lm(formula = SL ~ FemaleFK, data = MHWday15_30)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.3954 -0.2015 -0.1094 0.2004 0.6053
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.024509 0.183194 5.592 3.95e-08 ***
## FemaleFK 0.007113 0.034922 0.204 0.839
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2447 on 437 degrees of freedom
## (5 observations deleted due to missingness)
## Multiple R-squared: 9.492e-05, Adjusted R-squared: -0.002193
## F-statistic: 0.04149 on 1 and 437 DF, p-value: 0.8387lm = lm(SL ~ MaleSL, data=MHWday15_30)
summary(lm)##
## Call:
## lm(formula = SL ~ MaleSL, data = MHWday15_30)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.3486 -0.2099 -0.1026 0.2124 0.5882
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.23213 0.46074 4.845 1.76e-06 ***
## MaleSL -0.12155 0.04784 -2.541 0.0114 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2429 on 437 degrees of freedom
## (5 observations deleted due to missingness)
## Multiple R-squared: 0.01456, Adjusted R-squared: 0.01231
## F-statistic: 6.457 on 1 and 437 DF, p-value: 0.0114lm = lm(SL ~ MaleFK, data=MHWday15_30)
summary(lm)##
## Call:
## lm(formula = SL ~ MaleFK, data = MHWday15_30)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.37185 -0.20227 -0.09864 0.21132 0.59033
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.63872 0.13486 4.736 2.95e-06 ***
## MaleFK 0.08162 0.02592 3.148 0.00175 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.242 on 437 degrees of freedom
## (5 observations deleted due to missingness)
## Multiple R-squared: 0.02218, Adjusted R-squared: 0.01994
## F-statistic: 9.912 on 1 and 437 DF, p-value: 0.001754Strongest correlation with Female SL. Testing model fit.
lmer.SL15b = lmer(LogSL ~ Treatment * Day + (1|Clutch/Tank), data = MHWday15_30)
lmer.SL15c = lmer(LogSL ~ Treatment * Day + LogFemaleSL + (1|Clutch/Tank), data = MHWday15_30)
AIC(lmer.SL15b, lmer.SL15c)## df AIC
## lmer.SL15b 7 -958.6867
## lmer.SL15c 8 -962.5142Choosing c (with LogFemaleSL).
3. Final model
lmer.SL15 = lmer(LogSL ~ Treatment * Day + LogFemaleSL + (1|Clutch/Tank), data = MHWday15_30)4. Model performance
performance::check_model(lmer.SL15, check="homogeneity") 
performance::check_model(lmer.SL15, check="outliers") 
performance::check_model(lmer.SL15, check="qq", detrend = FALSE) 
performance::check_model(lmer.SL15, check="normality") 
performance::check_model(lmer.SL15, check="linearity") 
performance::check_model(lmer.SL15, check="pp_check")## Failed to compute posterior predictive checks with `re_formula=NULL`.
## Trying again with `re_formula=NA` now.
hist(residuals(lmer.SL15), col="darkgray")
shapiro.test(residuals(lmer.SL15))##
## Shapiro-Wilk normality test
##
## data: residuals(lmer.SL15)
## W = 0.97891, p-value = 5.345e-06qqline(resid(lmer.SL15))
outlierTest(lmer.SL15)## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
## rstudent unadjusted p-value Bonferroni p
## 257 3.667292 0.00027585 0.1211All outliers checked and removed if they were unusual for a given clutch.
5. Model outputs
Summary statistics
summary(lmer.SL15)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: LogSL ~ Treatment * Day + LogFemaleSL + (1 | Clutch/Tank)
## Data: MHWday15_30
##
## REML criterion at convergence: -978.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.4034 -0.5474 0.0168 0.6069 3.5624
##
## Random effects:
## Groups Name Variance Std.Dev.
## Tank:Clutch (Intercept) 0.004760 0.06899
## Clutch (Intercept) 0.001463 0.03825
## Residual 0.005164 0.07186
## Number of obs: 439, groups: Tank:Clutch, 24; Clutch, 6
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -4.57037 1.81595 4.01087 -2.517 0.0654 .
## TreatmentEarly exposure -0.00896 0.04101 14.97636 -0.218 0.8300
## Day30 0.40083 0.04102 14.98980 9.771 6.81e-08 ***
## LogFemaleSL 1.95658 0.80442 4.00907 2.432 0.0717 .
## TreatmentEarly exposure:Day30 -0.04056 0.05800 14.97620 -0.699 0.4951
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) TrtmEe Day30 LgFmSL
## TrtmntErlye -0.012
## Day30 -0.012 0.501
## LogFemaleSL -1.000 0.001 0.000
## TrtmnEe:D30 0.009 -0.707 -0.707 -0.001anova
Anova(lmer.SL15, type = 3)## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: LogSL
## Chisq Df Pr(>Chisq)
## (Intercept) 6.3342 1 0.01184 *
## Treatment 0.0477 1 0.82707
## Day 95.4693 1 < 2e-16 ***
## LogFemaleSL 5.9160 1 0.01500 *
## Treatment:Day 0.4889 1 0.48440
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Emmean
EmmSL.15_Log = (emmeans(lmer.SL15, ~ Treatment * Day) %>% as.data.frame)
# Back-transforming
EmmSL.15_30 <- EmmSL.15_Log %>%
mutate(
emmean = exp(emmean),
lower.CL = exp(lower.CL),
upper.CL = exp(upper.CL)
)
EmmSL.15_30## Treatment Day emmean SE df lower.CL upper.CL
## 1 Control 15 0.8588744 0.03297424 15.10866 0.8006179 0.9213700
## 2 Early exposure 15 0.8512136 0.03292342 15.01720 0.7935333 0.9130866
## 3 Control 30 1.2823491 0.03293410 15.03686 1.1954363 1.3755807
## 4 Early exposure 30 1.2203971 0.03296141 15.08607 1.1376397 1.3091748Pairwise
emm_day <- emmeans(lmer.SL15, ~ Day)## NOTE: Results may be misleading due to involvement in interactionspairs_day <- pairs(emm_day, adjust = "tukey") |> summary(infer = TRUE)
pairs_day## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Day15 - Day30 -0.381 0.029 15 -0.442 -0.319 -13.122 <.0001
##
## Results are averaged over the levels of: Treatment
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.956. Figure
fig.SL.15_30 <- ggplot(EmmSL.15_30,
aes(x = Day, y = emmean, fill = Treatment)) +
geom_errorbar(
aes(ymin = lower.CL, ymax = upper.CL),
position = position_dodge(width = 0.6),
width = 0.1, colour = "black", linewidth = 0.3
) +
geom_point(
aes(size = interaction(Treatment, Day)),
shape = 22, colour = "black", stroke = 0.3,
position = position_dodge(width = 0.6)
) +
scale_size_manual(
values = c(
"Control.15" = 6,
"Early exposure.15" = 6,
"Control.30" = 6,
"Early exposure.30" = 6
),
guide = "none"
) +
theme_classic() +
labs(
y = "Standard length (cm)"
) +
scale_fill_manual(values = c("Control" = "#0072B2",
"Early exposure" = "#F2C200"
)) +
scale_y_continuous(
limits = c(0.5, 1.5),
breaks = c(0.5, 1, 1.5),
labels = c(0.5, 1, 1.5)) +
theme(
text = element_text(family = "Helvetica"),
axis.title.x = element_blank(),
axis.title = element_text(size = 14, face = "bold", colour = "black"),
axis.text = element_text(size = 12, colour = "black"),
axis.text.x = element_blank(),
legend.position = "none"
)
fig.SL.15_30
Growth in weight (Day 15 and 30)
1. Random effect structures
lmer.W15a = lmer(LogW ~ Treatment * Day + (1|Tank) + (1|Clutch), data = MHWday15_30)
lmer.W15b = lmer(LogW ~ Treatment * Day + (1|Clutch/Tank) , data = MHWday15_30)
AIC(lmer.W15a, lmer.W15b)## df AIC
## lmer.W15a 7 51.58869
## lmer.W15b 7 50.73213Choosing b.
2. Covariate exploration
lm = lm(W ~ FemaleSL, data=MHWday15_30)
summary(lm)##
## Call:
## lm(formula = W ~ FemaleSL, data = MHWday15_30)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.04515 -0.02417 -0.01170 0.01435 0.11003
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.266890 0.059974 -4.450 1.09e-05 ***
## FemaleSL 0.033077 0.006265 5.279 2.05e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03317 on 436 degrees of freedom
## (6 observations deleted due to missingness)
## Multiple R-squared: 0.06009, Adjusted R-squared: 0.05793
## F-statistic: 27.87 on 1 and 436 DF, p-value: 2.049e-07lm = lm(W ~ FemaleFK, data=MHWday15_30)
summary(lm)##
## Call:
## lm(formula = W ~ FemaleFK, data = MHWday15_30)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.03916 -0.02550 -0.01460 0.01828 0.11194
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.064711 0.025614 2.526 0.0119 *
## FemaleFK -0.002882 0.004883 -0.590 0.5554
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0342 on 436 degrees of freedom
## (6 observations deleted due to missingness)
## Multiple R-squared: 0.0007982, Adjusted R-squared: -0.001494
## F-statistic: 0.3483 on 1 and 436 DF, p-value: 0.5554lm = lm(W ~ MaleSL, data=MHWday15_30)
summary(lm)##
## Call:
## lm(formula = W ~ MaleSL, data = MHWday15_30)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.03710 -0.02535 -0.01471 0.01874 0.11400
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.16706 0.06481 2.578 0.0103 *
## MaleSL -0.01220 0.00673 -1.812 0.0706 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03409 on 436 degrees of freedom
## (6 observations deleted due to missingness)
## Multiple R-squared: 0.007478, Adjusted R-squared: 0.005202
## F-statistic: 3.285 on 1 and 436 DF, p-value: 0.0706lm = lm(W ~ MaleFK, data=MHWday15_30)
summary(lm)##
## Call:
## lm(formula = W ~ MaleFK, data = MHWday15_30)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.03870 -0.02503 -0.01350 0.01852 0.11484
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.009677 0.018870 -0.513 0.60833
## MaleFK 0.011440 0.003627 3.154 0.00172 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03383 on 436 degrees of freedom
## (6 observations deleted due to missingness)
## Multiple R-squared: 0.02231, Adjusted R-squared: 0.02007
## F-statistic: 9.949 on 1 and 436 DF, p-value: 0.00172Strongest correlation with Female SL. Testing model fit.
lmer.W15b = lmer(LogW ~ Treatment * Day + (1|Clutch/Tank) , data = MHWday15_30)
lmer.W15c = lmer(LogW ~ Treatment * Day + LogFemaleSL + (1|Clutch/Tank) , data = MHWday15_30)
AIC(lmer.W15b, lmer.W15c)## df AIC
## lmer.W15b 7 50.73213
## lmer.W15c 8 45.34397Choosing c (with LogFemaleSL).
3. Final model
lmer.W15 = lmer(LogW ~ Treatment * Day + LogFemaleSL + (1|Clutch/Tank) , data = MHWday15_30)4. Model performance
performance::check_model(lmer.W15, check="homogeneity") 
performance::check_model(lmer.W15, check="outliers") 
performance::check_model(lmer.W15, check="qq", detrend = FALSE) 
performance::check_model(lmer.W15, check="normality") 
performance::check_model(lmer.W15, check="linearity") 
performance::check_model(lmer.W15, check="pp_check")## Failed to compute posterior predictive checks with `re_formula=NULL`.
## Trying again with `re_formula=NA` now.
hist(residuals(lmer.W15), col="darkgray")
shapiro.test(residuals(lmer.W15))##
## Shapiro-Wilk normality test
##
## data: residuals(lmer.W15)
## W = 0.98236, p-value = 3.587e-05qqline(resid(lmer.W15))
outlierTest(lmer.W15)## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
## rstudent unadjusted p-value Bonferroni p
## 143 -3.888549 0.00011683 0.051172#143 identified as an outlier but retained as it appears within the expected range for that clutch.
5. Model outputs
Summary statistics
summary(lmer.W15)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: LogW ~ Treatment * Day + LogFemaleSL + (1 | Clutch/Tank)
## Data: MHWday15_30
##
## REML criterion at convergence: 29.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.7615 -0.5614 -0.0094 0.6167 3.1889
##
## Random effects:
## Groups Name Variance Std.Dev.
## Tank:Clutch (Intercept) 0.06310 0.2512
## Clutch (Intercept) 0.03105 0.1762
## Residual 0.05192 0.2279
## Number of obs: 438, groups: Tank:Clutch, 24; Clutch, 6
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -18.25542 7.58553 4.00643 -2.407 0.0737
## TreatmentEarly exposure 0.08386 0.14831 15.01817 0.565 0.5801
## Day30 1.12433 0.14834 15.02843 7.579 1.64e-06
## LogFemaleSL 6.41820 3.36029 4.00511 1.910 0.1287
## TreatmentEarly exposure:Day30 -0.23208 0.20974 15.01502 -1.107 0.2859
##
## (Intercept) .
## TreatmentEarly exposure
## Day30 ***
## LogFemaleSL
## TreatmentEarly exposure:Day30
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) TrtmEe Day30 LgFmSL
## TrtmntErlye -0.010
## Day30 -0.010 0.501
## LogFemaleSL -1.000 0.001 0.000
## TrtmnEe:D30 0.008 -0.707 -0.707 -0.001anova
Anova(lmer.W15, type = 3)## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: LogW
## Chisq Df Pr(>Chisq)
## (Intercept) 5.7918 1 0.01610 *
## Treatment 0.3197 1 0.57180
## Day 57.4478 1 3.471e-14 ***
## LogFemaleSL 3.6482 1 0.05613 .
## Treatment:Day 1.2244 1 0.26850
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Emmean
EmmW.15_Log = emmeans(lmer.W15, ~ Treatment * Day) %>% as.data.frame()
# Back-transforming
EmmW.15_30 <- EmmW.15_Log %>%
mutate(
emmean = exp(emmean),
lower.CL = exp(lower.CL),
upper.CL = exp(upper.CL)
)
EmmW.15_30## Treatment Day emmean SE df lower.CL upper.CL
## 1 Control 15 0.02323872 0.1273055 12.94445 0.01764889 0.03059898
## 2 Early exposure 15 0.02527151 0.1271566 12.88449 0.01919639 0.03326924
## 3 Control 30 0.07153211 0.1271843 12.89579 0.05433429 0.09417336
## 4 Early exposure 30 0.06167778 0.1272555 12.92442 0.04684486 0.08120737Pairwise
emm_day <- emmeans(lmer.W15, ~ Day)## NOTE: Results may be misleading due to involvement in interactionspairs_day <- pairs(emm_day, adjust = "tukey") |> summary(infer = TRUE)
pairs_day## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Day15 - Day30 -1.01 0.105 15 -1.23 -0.785 -9.615 <.0001
##
## Results are averaged over the levels of: Treatment
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.956. Figure
fig.W.15_30 <- ggplot(EmmW.15_30,
aes(x = Day, y = emmean, fill = Treatment)) +
geom_errorbar(
aes(ymin = lower.CL, ymax = upper.CL),
position = position_dodge(width = 0.6),
width = 0.1, colour = "black", linewidth = 0.3
) +
geom_point(
aes(size = interaction(Treatment, Day)),
shape = 22, colour = "black", stroke = 0.3,
position = position_dodge(width = 0.6)
) +
scale_size_manual(
values = c(
"Control.15" = 6,
"Early exposure.15" = 6,
"Control.30" = 6,
"Early exposure.30" = 6
),
guide = "none"
) +
theme_classic() +
labs(
y = "Weight (g)"
) +
scale_fill_manual(values = c("Control" = "#0072B2",
"Early exposure" = "#F2C200"
)) +
scale_y_continuous(
limits = c(0, 0.1),
breaks = c(0, 0.05, 0.1),
labels = c(0, 0.05, 0.1)) +
theme(
text = element_text(family = "Helvetica"),
axis.title.x = element_blank(),
axis.title = element_text(size = 14, face = "bold", colour = "black"),
axis.text = element_text(size = 12, colour = "black"),
axis.text.x = element_blank(),
legend.position = "none"
)
fig.W.15_30
Growth in Condition (Day 15 and 30)
1. Random effect structure
MHWday15_30 <- MHWday15_30 %>%
group_by(Day) %>%
mutate(LogSL.centered = LogSL - mean(LogSL, na.rm = TRUE))lmer.Cond15a = lmer(LogW ~ Treatment * Day + LogSL.centered + (1|Tank) + (1|Clutch), data = MHWday15_30)
lmer.Cond15b = lmer(LogW ~ Treatment * Day + LogSL.centered + (1|Clutch/Tank), data = MHWday15_30)
AIC(lmer.Cond15a, lmer.Cond15b)## df AIC
## lmer.Cond15a 8 -779.3522
## lmer.Cond15b 8 -788.1354Choosing b.
2. Covariate exploration
Testing model fit with Female SL to follow the model for W above.
lmer.Cond15b = lmer(LogW ~ Treatment * Day + LogSL.centered + (1|Clutch/Tank), data = MHWday15_30)
lmer.Cond15c = lmer(LogW ~ Treatment * Day + LogSL.centered + LogFemaleSL + (1|Clutch/Tank), data = MHWday15_30)
AIC(lmer.Cond15b, lmer.Cond15c)## df AIC
## lmer.Cond15b 8 -788.1354
## lmer.Cond15c 9 -788.7881Choosing b (without covariate).
3. Final model
lmer.Cond15 = lmer(LogW ~ Treatment * Day + LogSL.centered + (1|Clutch/Tank), data = MHWday15_30)4. Model performance
performance::check_model(lmer.Cond15, check="homogeneity") 
performance::check_model(lmer.Cond15, check="outliers") 
performance::check_model(lmer.Cond15, check="qq", detrend = FALSE) 
performance::check_model(lmer.Cond15, check="normality") 
performance::check_model(lmer.Cond15, check="linearity") 
performance::check_model(lmer.Cond15, check="pp_check")## Failed to compute posterior predictive checks with `re_formula=NULL`.
## Trying again with `re_formula=NA` now.
hist(residuals(lmer.Cond15), col="darkgray")
shapiro.test(residuals(lmer.Cond15))##
## Shapiro-Wilk normality test
##
## data: residuals(lmer.Cond15)
## W = 0.99544, p-value = 0.2303qqline(resid(lmer.Cond15))
outlierTest(lmer.Cond15)## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
## rstudent unadjusted p-value Bonferroni p
## 339 -3.826096 0.00014953 0.065496#339 identified as an outlier but retained as it appears within the expected range of the clutch.
5. Model output
Summary statistics
summary(lmer.Cond15)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: LogW ~ Treatment * Day + LogSL.centered + (1 | Clutch/Tank)
## Data: MHWday15_30
##
## REML criterion at convergence: -804.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.7057 -0.6139 0.0062 0.5649 2.9462
##
## Random effects:
## Groups Name Variance Std.Dev.
## Tank:Clutch (Intercept) 0.016541 0.12861
## Clutch (Intercept) 0.002562 0.05062
## Residual 0.007312 0.08551
## Number of obs: 438, groups: Tank:Clutch, 24; Clutch, 6
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -3.78711 0.05705 18.91185 -66.382 < 2e-16
## TreatmentEarly exposure 0.11381 0.07517 15.00680 1.514 0.151
## Day30 1.09000 0.07517 15.01288 14.500 3.09e-10
## LogSL.centered 2.95780 0.05780 428.00498 51.170 < 2e-16
## TreatmentEarly exposure:Day30 -0.11549 0.10632 15.01669 -1.086 0.295
##
## (Intercept) ***
## TreatmentEarly exposure
## Day30 ***
## LogSL.centered ***
## TreatmentEarly exposure:Day30
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) TrtmEe Day30 LgSL.c
## TrtmntErlye -0.659
## Day30 -0.659 0.500
## LogSL.cntrd -0.006 0.008 -0.009
## TrtmnEe:D30 0.466 -0.707 -0.707 0.021anova
Anova(lmer.Cond15, type = 3)## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: LogW
## Chisq Df Pr(>Chisq)
## (Intercept) 4406.5041 1 <2e-16 ***
## Treatment 2.2924 1 0.1300
## Day 210.2369 1 <2e-16 ***
## LogSL.centered 2618.4124 1 <2e-16 ***
## Treatment:Day 1.1799 1 0.2774
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Emmean
EmmCond.15_Log = (emmeans(lmer.Cond15, ~ Treatment * Day) %>% as.data.frame)
# Back-transforming
EmmCond.15_30 <- EmmCond.15_Log %>%
mutate(
emmean = exp(emmean),
lower.CL = exp(lower.CL),
upper.CL = exp(upper.CL)
)
EmmCond.15_30## Treatment Day emmean SE df lower.CL upper.CL
## 1 Control 15 0.02269324 0.05705064 18.98494 0.02013888 0.02557159
## 2 Early exposure 15 0.02542861 0.05700263 18.92156 0.02256802 0.02865181
## 3 Control 30 0.06749562 0.05701975 18.94402 0.05990112 0.07605299
## 4 Early exposure 30 0.06738237 0.05706555 19.00242 0.05979636 0.07593076Pairwise
Day
emm_day <- emmeans(lmer.Cond15, ~ Day)## NOTE: Results may be misleading due to involvement in interactionspairs_day <- pairs(emm_day, adjust = "tukey") |> summary(infer = TRUE)
pairs_day## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Day15 - Day30 -1.03 0.0532 15 -1.15 -0.919 -19.421 <.0001
##
## Results are averaged over the levels of: Treatment
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.956. Figure
fig.Cond.15_30 <- ggplot(EmmCond.15_30,
aes(x = Day, y = emmean, fill = Treatment)) +
geom_errorbar(
aes(ymin = lower.CL, ymax = upper.CL),
position = position_dodge(width = 0.6),
width = 0.1, colour = "black", linewidth = 0.3
) +
geom_point(
aes(size = interaction(Treatment, Day)),
shape = 22, colour = "black", stroke = 0.3,
position = position_dodge(width = 0.6)
) +
scale_size_manual(
values = c(
"Control.15" = 4.5,
"Early exposure.15" = 4.5,
"Control.30" = 6,
"Early exposure.30" = 6
),
guide = "none"
) +
theme_classic() +
labs(
y = "Condition (g)"
) +
scale_fill_manual(values = c("Control" = "#0072B2",
"Early exposure" = "#F2C200"
)) +
scale_y_continuous(
limits = c(0.02, 0.08),
breaks = c(0.02, 0.04, 0.06, 0.08),
labels = c(0.02, 0.04, 0.06, 0.08)) +
theme(
text = element_text(family = "Helvetica"),
axis.title.x = element_blank(),
axis.title = element_text(size = 14, face = "bold", colour = "black"),
axis.text = element_text(size = 12, colour = "black"),
legend.position = "none"
)
fig.Cond.15_30
Growth in SL (Day 45 and 60)
1. Random effect structures
lmer.SL45a = lmer(SL ~ Treatment * Day + (1|Tank) + (1|Clutch), data = MHWday45_60)
lmer.SL45b = lmer(SL ~ Treatment * Day + (1|Clutch/Tank), data = MHWday45_60)
AIC(lmer.SL45a, lmer.SL45b)## df AIC
## lmer.SL45a 9 -412.6339
## lmer.SL45b 9 -415.3363Choosing b (without covariate).
2. Covariate exploration
lm = lm(SL ~ FemaleSL, data=MHWday45_60)
summary(lm)##
## Call:
## lm(formula = SL ~ FemaleSL, data = MHWday45_60)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.69332 -0.18732 -0.01752 0.20068 0.65450
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.92359 0.40934 2.256 0.02442 *
## FemaleSL 0.11211 0.04279 2.620 0.00902 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2608 on 591 degrees of freedom
## (9 observations deleted due to missingness)
## Multiple R-squared: 0.01148, Adjusted R-squared: 0.009809
## F-statistic: 6.864 on 1 and 591 DF, p-value: 0.00902lm = lm(SL ~ FemaleFK, data=MHWday45_60)
summary(lm)##
## Call:
## lm(formula = SL ~ FemaleFK, data = MHWday45_60)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.6688 -0.1860 -0.0198 0.2014 0.6246
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.77668 0.17192 10.335 <2e-16 ***
## FemaleFK 0.04206 0.03296 1.276 0.202
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.262 on 591 degrees of freedom
## (9 observations deleted due to missingness)
## Multiple R-squared: 0.002749, Adjusted R-squared: 0.001061
## F-statistic: 1.629 on 1 and 591 DF, p-value: 0.2023lm = lm(SL ~ MaleSL, data=MHWday45_60)
summary(lm)##
## Call:
## lm(formula = SL ~ MaleSL, data = MHWday45_60)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.69531 -0.18992 -0.01028 0.19869 0.62072
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.23371 0.42133 0.555 0.579
## MaleSL 0.18272 0.04368 4.183 3.31e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2585 on 591 degrees of freedom
## (9 observations deleted due to missingness)
## Multiple R-squared: 0.02876, Adjusted R-squared: 0.02711
## F-statistic: 17.5 on 1 and 591 DF, p-value: 3.31e-05lm = lm(SL ~ MaleFK, data=MHWday45_60)
summary(lm)##
## Call:
## lm(formula = SL ~ MaleFK, data = MHWday45_60)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.69030 -0.19254 -0.00554 0.19529 0.63001
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.32821 0.13195 17.645 <2e-16 ***
## MaleFK -0.06472 0.02559 -2.529 0.0117 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2609 on 591 degrees of freedom
## (9 observations deleted due to missingness)
## Multiple R-squared: 0.0107, Adjusted R-squared: 0.009029
## F-statistic: 6.394 on 1 and 591 DF, p-value: 0.01171Correlation strongest with Male SL Testing model fit.
lmer.SL45b = lmer(SL ~ Treatment * Day + (1|Clutch/Tank), data = MHWday45_60)
lmer.SL45c = lmer(SL ~ Treatment * Day + MaleSL + (1|Clutch/Tank), data = MHWday45_60)## boundary (singular) fit: see help('isSingular')AIC(lmer.SL45b, lmer.SL45c)## df AIC
## lmer.SL45b 9 -415.3363
## lmer.SL45c 10 -413.7352Choosing b (without covariate).
3. Final model
lmer.SL45 = lmer(SL ~ Treatment * Day + (1|Clutch/Tank), data = MHWday45_60)4. Performance check
performance::check_model(lmer.SL45, check="homogeneity") 
performance::check_model(lmer.SL45, check="outliers") 
performance::check_model(lmer.SL45, check="qq", detrend = FALSE) 
performance::check_model(lmer.SL45, check="normality") 
performance::check_model(lmer.SL45, check="linearity") 
performance::check_model(lmer.SL45, check="pp_check")## Failed to compute posterior predictive checks with `re_formula=NULL`.
## Trying again with `re_formula=NA` now.
hist(residuals(lmer.SL45), col="darkgray")
qqline(resid(lmer.SL45))
outlierTest(lmer.SL45)## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
## rstudent unadjusted p-value Bonferroni p
## 80 3.547752 0.0004198 0.24894#80 identified as an outlier but retained as it appears within expected range for that clutch.
5. Outputs
Summary
summary(lmer.SL45)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: SL ~ Treatment * Day + (1 | Clutch/Tank)
## Data: MHWday45_60
##
## REML criterion at convergence: -433.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1231 -0.6408 0.0319 0.6786 3.4604
##
## Random effects:
## Groups Name Variance Std.Dev.
## Tank:Clutch (Intercept) 0.011738 0.10834
## Clutch (Intercept) 0.001234 0.03513
## Residual 0.023905 0.15461
## Number of obs: 593, groups: Tank:Clutch, 34; Clutch, 6
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.77997 0.05345 27.32428 33.301 < 2e-16 ***
## TreatmentEarly exposure 0.04458 0.06947 22.93183 0.642 0.527
## TreatmentLate exposure 0.06753 0.06943 22.86459 0.973 0.341
## Day60 0.41974 0.07318 23.99935 5.736 6.55e-06 ***
## TreatmentEarly exposure:Day60 -0.03385 0.09866 23.34803 -0.343 0.735
## TreatmentLate exposure:Day60 -0.09779 0.09859 23.29326 -0.992 0.331
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) TrtmEe TrtmLe Day60 TEe:D6
## TrtmntErlye -0.714
## TrtmntLtexp -0.715 0.550
## Day60 -0.679 0.523 0.523
## TrtmnEe:D60 0.504 -0.705 -0.388 -0.742
## TrtmnLe:D60 0.504 -0.388 -0.705 -0.742 0.550anova
Anova(lmer.SL45, type = 3)## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: SL
## Chisq Df Pr(>Chisq)
## (Intercept) 1108.9536 1 < 2.2e-16 ***
## Treatment 0.9625 2 0.6180
## Day 32.8997 1 9.704e-09 ***
## Treatment:Day 1.0429 2 0.5936
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Emmeans
EmmSL.45_60 = (emmeans(lmer.SL45, ~ Treatment * Day) %>% as.data.frame)
EmmSL.45_60## Treatment Day emmean SE df lower.CL upper.CL
## Control 45 1.779972 0.05391655 27.43 1.669425 1.890519
## Early exposure 45 1.824557 0.04879200 26.71 1.724394 1.924720
## Late exposure 45 1.847500 0.04872323 26.57 1.747452 1.947548
## Control 60 2.199708 0.05428961 28.19 2.088535 2.310881
## Early exposure 60 2.210445 0.04911352 27.40 2.109741 2.311149
## Late exposure 60 2.169442 0.04902208 27.14 2.068882 2.270002
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95Pairwise
Day
emm_day <- emmeans(lmer.SL45, ~ Day)## NOTE: Results may be misleading due to involvement in interactionspairs_day <- pairs(emm_day, adjust = "tukey") |> summary(infer = TRUE)
pairs_day## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Day45 - Day60 -0.376 0.0398 23.5 -0.458 -0.294 -9.454 <.0001
##
## Results are averaged over the levels of: Treatment
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.956. Figure
fig.SL.45_60 <- ggplot(EmmSL.45_60,
aes(x = Day, y = emmean, fill = Treatment)) +
geom_errorbar(
aes(ymin = lower.CL, ymax = upper.CL),
position = position_dodge(width = 0.6),
width = 0.1, colour = "black", linewidth = 0.3
) +
geom_point(
aes(size = interaction(Treatment, Day)),
shape = 22, colour = "black", stroke = 0.3,
position = position_dodge(width = 0.6)
) +
scale_size_manual(
values = c(
"Control.45" = 6,
"Early exposure.45" = 6,
"Late exposure.45" = 6,
"Control.60" = 6,
"Early exposure.60" = 6,
"Late exposure.60" = 6
),
guide = "none"
) +
theme_classic() +
labs(
y = "Standard length (cm)"
) +
scale_fill_manual(values = c("Control" = "#0072B2",
"Early exposure" = "#F2C200",
"Late exposure" = "#D55E00")) +
scale_y_continuous(
limits = c(1.5, 2.5),
breaks = c(1.5, 2, 2.5),
labels = c(1.5, 2, 2.5)) +
theme(
text = element_text(family = "Helvetica"),
axis.title.x = element_blank(),
axis.title = element_text(size = 14, face = "bold", colour = "black"),
axis.text = element_text(size = 12, colour = "black"),
axis.text.x = element_blank(),
legend.position = "none"
)
fig.SL.45_60
# Grwoth in weight (Day45 and 60) ## 1. Random effect structures
lmer.W45a = lmer(W ~ Treatment * Day + (1|Tank) + (1|Clutch), data = MHWday45_60)
lmer.W45b = lmer(W ~ Treatment * Day + (1|Clutch/Tank), data = MHWday45_60)
AIC(lmer.W45a, lmer.W45b)## df AIC
## lmer.W45a 9 -1287.555
## lmer.W45b 9 -1287.433Choosing b.
2. Covaraite exploration
lm = lm(W ~ FemaleSL, data=MHWday45_60)
summary(lm)##
## Call:
## lm(formula = W ~ FemaleSL, data = MHWday45_60)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.19801 -0.08580 -0.02011 0.07576 0.43739
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.13159 0.17933 0.734 0.463
## FemaleSL 0.01783 0.01875 0.951 0.342
##
## Residual standard error: 0.1143 on 592 degrees of freedom
## (8 observations deleted due to missingness)
## Multiple R-squared: 0.001526, Adjusted R-squared: -0.0001604
## F-statistic: 0.9049 on 1 and 592 DF, p-value: 0.3419lm = lm(W ~ FemaleFK, data=MHWday45_60)
summary(lm)##
## Call:
## lm(formula = W ~ FemaleFK, data = MHWday45_60)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.20780 -0.08380 -0.01960 0.07286 0.44326
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.23607 0.07493 3.151 0.00171 **
## FemaleFK 0.01269 0.01437 0.883 0.37749
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1143 on 592 degrees of freedom
## (8 observations deleted due to missingness)
## Multiple R-squared: 0.001316, Adjusted R-squared: -0.0003711
## F-statistic: 0.78 on 1 and 592 DF, p-value: 0.3775lm = lm(W ~ MaleSL, data=MHWday45_60)
summary(lm)##
## Call:
## lm(formula = W ~ MaleSL, data = MHWday45_60)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.19595 -0.08816 -0.02220 0.07690 0.43475
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.30264 0.18470 -1.639 0.10184
## MaleSL 0.06272 0.01915 3.275 0.00112 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1133 on 592 degrees of freedom
## (8 observations deleted due to missingness)
## Multiple R-squared: 0.0178, Adjusted R-squared: 0.01614
## F-statistic: 10.73 on 1 and 592 DF, p-value: 0.001117lm = lm(W ~ MaleFK, data=MHWday45_60)
summary(lm)##
## Call:
## lm(formula = W ~ MaleFK, data = MHWday45_60)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.19489 -0.08481 -0.01922 0.07493 0.43581
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.44485 0.05753 7.733 4.54e-14 ***
## MaleFK -0.02778 0.01116 -2.489 0.0131 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1138 on 592 degrees of freedom
## (8 observations deleted due to missingness)
## Multiple R-squared: 0.01036, Adjusted R-squared: 0.008688
## F-statistic: 6.197 on 1 and 592 DF, p-value: 0.01307Correlation with Male SL. Testing model fit.
lmer.W45b = lmer(W ~ Treatment * Day + (1|Clutch/Tank), data = MHWday45_60)
lmer.W45c = lmer(W ~ Treatment * Day + MaleSL + (1|Clutch/Tank), data = MHWday45_60)
AIC(lmer.W45b, lmer.W45c)## df AIC
## lmer.W45b 9 -1287.433
## lmer.W45c 10 -1282.685Choosing b (without covariate).
3. Final model
lmer.W45 = lmer(W ~ Treatment * Day + (1|Clutch/Tank), data = MHWday45_60)4. Performance check
performance::check_model(lmer.W45, check="homogeneity") 
performance::check_model(lmer.W45, check="outliers") 
performance::check_model(lmer.W45, check="qq", detrend = FALSE) 
performance::check_model(lmer.W45, check="normality") 
performance::check_model(lmer.W45, check="linearity") 
performance::check_model(lmer.W45, check="pp_check")## Failed to compute posterior predictive checks with `re_formula=NULL`.
## Trying again with `re_formula=NA` now.
hist(residuals(lmer.W45), col="darkgray")
shapiro.test(residuals(lmer.W45))##
## Shapiro-Wilk normality test
##
## data: residuals(lmer.W45)
## W = 0.98056, p-value = 4.213e-07qqline(resid(lmer.W45))
outlierTest(lmer.W45)## rstudent unadjusted p-value Bonferroni p
## 79 5.434820 8.0666e-08 4.7915e-05
## 80 4.605901 5.0445e-06 2.9964e-03Outlier appears within the range of clutches. Kept.
5. Outputs
Summary statistic
summary(lmer.W45)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: W ~ Treatment * Day + (1 | Clutch/Tank)
## Data: MHWday45_60
##
## REML criterion at convergence: -1305.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.8026 -0.6200 -0.0684 0.6034 5.3113
##
## Random effects:
## Groups Name Variance Std.Dev.
## Tank:Clutch (Intercept) 0.0018655 0.04319
## Clutch (Intercept) 0.0001471 0.01213
## Residual 0.0055224 0.07431
## Number of obs: 594, groups: Tank:Clutch, 34; Clutch, 6
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.149e-01 2.156e-02 2.663e+01 9.968 1.76e-10
## TreatmentEarly exposure 1.545e-02 2.832e-02 2.202e+01 0.546 0.591
## TreatmentLate exposure 2.114e-02 2.829e-02 2.193e+01 0.747 0.463
## Day60 1.634e-01 2.986e-02 2.326e+01 5.471 1.41e-05
## TreatmentEarly exposure:Day60 6.135e-05 4.026e-02 2.252e+01 0.002 0.999
## TreatmentLate exposure:Day60 -1.528e-02 4.022e-02 2.245e+01 -0.380 0.708
##
## (Intercept) ***
## TreatmentEarly exposure
## TreatmentLate exposure
## Day60 ***
## TreatmentEarly exposure:Day60
## TreatmentLate exposure:Day60
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) TrtmEe TrtmLe Day60 TEe:D6
## TrtmntErlye -0.721
## TrtmntLtexp -0.722 0.550
## Day60 -0.685 0.522 0.522
## TrtmnEe:D60 0.508 -0.704 -0.387 -0.742
## TrtmnLe:D60 0.509 -0.387 -0.704 -0.742 0.551anova
Anova(lmer.W45, type = 3)## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: W
## Chisq Df Pr(>Chisq)
## (Intercept) 99.3532 1 < 2.2e-16 ***
## Treatment 0.5842 2 0.7467
## Day 29.9349 1 4.468e-08 ***
## Treatment:Day 0.2080 2 0.9012
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Emmeans
EmmW.45_60 = (emmeans(lmer.W45, ~ Treatment * Day) %>% as.data.frame)
EmmW.45_60## Treatment Day emmean SE df lower.CL upper.CL
## Control 45 0.2149288 0.02177094 27.63 0.1703059 0.2595518
## Early exposure 45 0.2303810 0.01964464 27.00 0.1900733 0.2706886
## Late exposure 45 0.2360643 0.01960534 26.79 0.1958226 0.2763061
## Control 60 0.3783079 0.02198242 28.71 0.3333290 0.4232868
## Early exposure 60 0.3938214 0.01981646 27.89 0.3532222 0.4344206
## Late exposure 60 0.3841627 0.01977142 27.55 0.3436329 0.4246925
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95Pairwise
Day
emm_day <- emmeans(lmer.W45, ~ Day)## NOTE: Results may be misleading due to involvement in interactionspairs_day <- pairs(emm_day, adjust = "tukey") |> summary(infer = TRUE)
pairs_day## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Day45 - Day60 -0.158 0.0162 23.5 -0.192 -0.125 -9.756 <.0001
##
## Results are averaged over the levels of: Treatment
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.956. Figure
fig.W.45_60 <- ggplot(EmmW.45_60,
aes(x = Day, y = emmean, fill = Treatment)) +
geom_errorbar(
aes(ymin = lower.CL, ymax = upper.CL),
position = position_dodge(width = 0.6),
width = 0.1, colour = "black", linewidth = 0.3
) +
geom_point(
aes(size = interaction(Treatment, Day)),
shape = 22, colour = "black", stroke = 0.3,
position = position_dodge(width = 0.6)
) +
scale_size_manual(
values = c(
"Control.45" = 6,
"Early exposure.45" = 6,
"Late exposure.45" = 6,
"Control.60" = 6,
"Early exposure.60" = 6,
"Late exposure.60" = 6
),
guide = "none"
) +
theme_classic() +
labs(
y = "Weight (g)"
) +
scale_fill_manual(values = c("Control" = "#0072B2",
"Early exposure" = "#F2C200",
"Late exposure" = "#D55E00")) +
scale_y_continuous(
limits = c(0.1, 0.5),
breaks = c(0.1, 0.3, 0.5),
labels = c(0.1, 0.3, 0.5)) +
theme(
text = element_text(family = "Helvetica"),
axis.title.x = element_blank(),
axis.title = element_text(size = 14, face = "bold", colour = "black"),
axis.text = element_text(size = 12, colour = "black"),
axis.text.x = element_blank(),
legend.position = "none"
)
fig.W.45_60
Growth in Condition (Day 45 and 60)
1. Random effect structure
Exploring whether SL and W data should be logged.
ggplot(MHWday45_60, aes(x = SL, y = W)) +
geom_point() +
geom_smooth(method = "lm", col = "red") +
ggtitle("Scatter Plot with Best Fit Line") +
xlab("SL") +
ylab("W")## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 9 rows containing non-finite outside the scale range
## (`stat_smooth()`).## Warning: Removed 9 rows containing missing values or values outside the scale range
## (`geom_point()`).
ggplot(MHWday45_60, aes(x = SL, y = W)) +
geom_point() +
geom_smooth(method = "lm", formula = y ~ poly(x, 3), col = "blue") +
ggtitle("Scatter Plot with Cubic Fit Line") +
xlab("SL") +
ylab("W")## Warning: Removed 9 rows containing non-finite outside the scale range (`stat_smooth()`).
## Removed 9 rows containing missing values or values outside the scale range
## (`geom_point()`).
ggplot(MHWday45_60, aes(x = LogSL, y = LogW)) +
geom_point() +
geom_smooth(method = "lm", col = "red") +
ggtitle("Scatter Plot with Best Fit Line") +
xlab("LogSL") +
ylab("LogW")## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 9 rows containing non-finite outside the scale range (`stat_smooth()`).
## Removed 9 rows containing missing values or values outside the scale range
## (`geom_point()`).
ggplot(MHWday45_60, aes(x = LogSL, y = LogW)) +
geom_point() +
geom_smooth(method = "lm", formula = y ~ poly(x, 3), col = "blue") +
ggtitle("Scatter Plot with Cubic Fit Line") +
xlab("LogSL") +
ylab("LogW")## Warning: Removed 9 rows containing non-finite outside the scale range (`stat_smooth()`).
## Removed 9 rows containing missing values or values outside the scale range
## (`geom_point()`).
Logging data to meet linearity.
MHWday45_60 <- MHWday45_60 %>%
group_by(Day) %>%
mutate(LogSL.centered = LogSL - mean(LogSL, na.rm = TRUE))lmer.Cond45a = lmer(LogW ~ Treatment * Day + LogSL.centered + (1|Tank) + (1|Clutch), data = MHWday45_60)
lmer.Cond45b = lmer(LogW ~ Treatment * Day + LogSL.centered + (1|Clutch/Tank), data = MHWday45_60)
AIC(lmer.Cond45a, lmer.Cond45b)## df AIC
## lmer.Cond45a 10 -1440.123
## lmer.Cond45b 10 -1461.438Choosing b.
2. Covariate exploration
Following weight, testing model fit including Male SL.
lmer.Cond45b = lmer(LogW ~ Treatment * Day + LogSL.centered + (1|Clutch/Tank), data = MHWday45_60)
lmer.Cond45c = lmer(LogW ~ Treatment * Day + LogSL.centered + LogMaleSL + (1|Clutch/Tank), data = MHWday45_60)
AIC(lmer.Cond45b, lmer.Cond45c)## df AIC
## lmer.Cond45b 10 -1461.438
## lmer.Cond45c 11 -1461.202Choosing b (without covariate).
3. Final model
lmer.Cond45 = lmer(LogW ~ Treatment * Day + LogSL.centered + (1|Clutch/Tank), data = MHWday45_60)4. Model performance
performance::check_model(lmer.Cond45, check="homogeneity") 
performance::check_model(lmer.Cond45, check="outliers") 
performance::check_model(lmer.Cond45, check="qq", detrend = FALSE) 
performance::check_model(lmer.Cond45, check="normality") 
performance::check_model(lmer.Cond45, check="linearity") 
performance::check_model(lmer.Cond45, check="pp_check")
hist(residuals(lmer.Cond45), col="darkgray")
shapiro.test(residuals(lmer.Cond45))##
## Shapiro-Wilk normality test
##
## data: residuals(lmer.Cond45)
## W = 0.96814, p-value = 4.745e-10qqline(resid(lmer.Cond45))
outlierTest(lmer.Cond45)## rstudent unadjusted p-value Bonferroni p
## 516 6.522674 1.5014e-10 8.9033e-08
## 528 4.746657 2.6061e-06 1.5454e-03Both outlier appear to be within the expected range of a clutch. Kept.
5. Outputs
Summary statistics
summary(lmer.Cond45)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: LogW ~ Treatment * Day + LogSL.centered + (1 | Clutch/Tank)
## Data: MHWday45_60
##
## REML criterion at convergence: -1481.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.4917 -0.6294 -0.0324 0.5969 6.3224
##
## Random effects:
## Groups Name Variance Std.Dev.
## Tank:Clutch (Intercept) 0.0029673 0.05447
## Clutch (Intercept) 0.0005685 0.02384
## Residual 0.0038979 0.06243
## Number of obs: 593, groups: Tank:Clutch, 34; Clutch, 6
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -1.516132 0.027272 26.304910 -55.594 < 2e-16
## TreatmentEarly exposure -0.008564 0.034341 23.413608 -0.249 0.805
## TreatmentLate exposure -0.002685 0.034341 23.407644 -0.078 0.938
## Day60 0.494959 0.036177 24.197670 13.682 6.94e-13
## LogSL.centered 2.989828 0.032945 582.559337 90.753 < 2e-16
## TreatmentEarly exposure:Day60 0.022549 0.048718 23.697997 0.463 0.648
## TreatmentLate exposure:Day60 0.054379 0.048724 23.710235 1.116 0.276
##
## (Intercept) ***
## TreatmentEarly exposure
## TreatmentLate exposure
## Day60 ***
## LogSL.centered ***
## TreatmentEarly exposure:Day60
## TreatmentLate exposure:Day60
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) TrtmEe TrtmLe Day60 LgSL.c TEe:D6
## TrtmntErlye -0.693
## TrtmntLtexp -0.693 0.551
## Day60 -0.660 0.524 0.525
## LogSL.cntrd 0.027 -0.022 -0.038 -0.025
## TrtmnEe:D60 0.490 -0.706 -0.389 -0.742 0.012
## TrtmnLe:D60 0.491 -0.390 -0.706 -0.743 0.035 0.551anova
Anova(lmer.Cond45, type = 3)## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: LogW
## Chisq Df Pr(>Chisq)
## (Intercept) 3090.6607 1 <2e-16 ***
## Treatment 0.0672 2 0.9670
## Day 187.1912 1 <2e-16 ***
## LogSL.centered 8236.1500 1 <2e-16 ***
## Treatment:Day 1.2790 2 0.5275
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Emmeans
EmmCond.45_60.Log = (emmeans(lmer.Cond45, ~ Treatment * Day) %>% as.data.frame)
EmmCond.45_60 <- EmmCond.45_60.Log %>%
mutate(
emmean = exp(emmean),
lower.CL = exp(lower.CL),
upper.CL = exp(upper.CL)
)
EmmCond.45_60## Treatment Day emmean SE df lower.CL upper.CL
## 1 Control 45 0.2195595 0.02745030 26.40380 0.2075226 0.2322945
## 2 Early exposure 45 0.2176872 0.02499873 25.05850 0.2067643 0.2291872
## 3 Late exposure 45 0.2189708 0.02498296 24.99540 0.2079889 0.2305326
## 4 Control 60 0.3601723 0.02756187 26.83081 0.3403633 0.3811342
## 5 Early exposure 60 0.3652448 0.02510444 25.47442 0.3468572 0.3846073
## 6 Late exposure 60 0.3792809 0.02507590 25.34261 0.3602030 0.3993693Pairwise
Day
emm_day <- emmeans(lmer.Cond45, ~ Day)## NOTE: Results may be misleading due to involvement in interactionspairs_day <- pairs(emm_day, adjust = "tukey") |> summary(infer = TRUE)
pairs_day## contrast estimate SE df lower.CL upper.CL t.ratio p.value
## Day45 - Day60 -0.521 0.0196 23.3 -0.561 -0.48 -26.548 <.0001
##
## Results are averaged over the levels of: Treatment
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.956. Figure
fig.Cond.45_60 <- ggplot(EmmCond.45_60,
aes(x = Day, y = emmean, fill = Treatment)) +
geom_errorbar(
aes(ymin = lower.CL, ymax = upper.CL),
position = position_dodge(width = 0.6),
width = 0.1, colour = "black", linewidth = 0.3
) +
geom_point(
aes(size = interaction(Treatment, Day)),
shape = 22, colour = "black", stroke = 0.3,
position = position_dodge(width = 0.6)
) +
scale_size_manual(
values = c(
"Control.45" = 6,
"Early exposure.45" = 6,
"Late exposure.45" = 6,
"Control.60" = 6,
"Early exposure.60" = 6,
"Late exposure.60" = 6
),
guide = "none"
) +
theme_classic() +
labs(
y = "Condition (g)"
) +
scale_fill_manual(values = c("Control" = "#0072B2",
"Early exposure" = "#F2C200",
"Late exposure" = "#D55E00")) +
scale_y_continuous(
limits = c(0.175, 0.4),
breaks = c(0.2, 0.3, 0.4),
labels = c(0.2, 0.3, 0.4)) +
theme(
text = element_text(family = "Helvetica"),
axis.title.x = element_blank(),
axis.title = element_text(size = 14, face = "bold", colour = "black"),
axis.text = element_text(size = 12, colour = "black"),
legend.position = "none"
)
fig.Cond.45_60