Code - Response of early life-stages of forest-forming seaweeds from warm-edge and central populations to marine heatwaves
2025-08-12
This code include the steps to obtain figures and results (statistical analyses) from the manuscript “Response of early life-stages of forest-forming seaweeds from warm-edge and central populations to marine heatwave”. It should be reproducible if the all data is stored in a folder called “Data_warming”.
Clear the environment
rm(list = ls())Load packages
Before loading the packages install them using install()
library(summarytools)
library(tidyr)
library(dplyr)
library(patchwork)
library(ggplot2)
library(lme4)
library(emmeans)
library(car)
library(outliers)
#library(vegan)
library(skimr)
library(gt)
library(RVAideMemoire)
library(pander)Section 1: Import, examine and manipulate data
In this section data frames were inspected to check variable names and types. In some situations the dataset was slightly modified for further manipulation, and summary statistics were calculated.
1a) Import the data
Count_germlings <- read.csv("Data_warming/Zygotes_counts.csv")
Length_germlings <- read.csv("Data_warming/Zygotes_length.csv")
Juveniles <- read.csv("Data_warming/Measurements_juveniles1_analyses.csv")
Tanks_metadata <- read.csv("Data_warming/Tanks_metadata.csv")
Acetone <- read.csv("Data_warming/Acetone.csv")
DMSO <- read.csv("Data_warming/DMSO.csv")
Pigments <- read.csv("Data_warming/Pigments.csv")
CN_results <- read.csv("Data_warming/CN_results.csv")1b) Check data structure and summarise dataset
1. Germling counts
- Check data structure
Count_germlings %>%
skimr::skim() %>%
gt::gt()| skim_type | skim_variable | n_missing | complete_rate | character.min | character.max | character.empty | character.n_unique | character.whitespace | numeric.mean | numeric.sd | numeric.p0 | numeric.p25 | numeric.p50 | numeric.p75 | numeric.p100 | numeric.hist |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| character | Treatment_Code | 0 | 1 | 2 | 2 | 0 | 3 | 0 | NA | NA | NA | NA | NA | NA | NA | NA |
| character | Treatment | 0 | 1 | 4 | 7 | 0 | 3 | 0 | NA | NA | NA | NA | NA | NA | NA | NA |
| character | Plate | 0 | 1 | 1 | 1 | 0 | 2 | 0 | NA | NA | NA | NA | NA | NA | NA | NA |
| character | Site | 0 | 1 | 7 | 12 | 0 | 6 | 0 | NA | NA | NA | NA | NA | NA | NA | NA |
| character | Date | 0 | 1 | 9 | 10 | 0 | 8 | 0 | NA | NA | NA | NA | NA | NA | NA | NA |
| numeric | Tank | 0 | 1 | NA | NA | NA | NA | NA | 3.500000 | 1.708320 | 1 | 2.00 | 3.5 | 5.00 | 6 | ▇▃▃▃▃ |
| numeric | Number | 0 | 1 | NA | NA | NA | NA | NA | 3.500000 | 1.708320 | 1 | 2.00 | 3.5 | 5.00 | 6 | ▇▃▃▃▃ |
| numeric | Day | 0 | 1 | NA | NA | NA | NA | NA | 21.875000 | 8.697882 | 9 | 15.25 | 21.5 | 27.75 | 37 | ▇▇▃▇▃ |
| numeric | Count_Alive | 0 | 1 | NA | NA | NA | NA | NA | 1.487269 | 2.008626 | 0 | 0.00 | 1.0 | 2.00 | 11 | ▇▂▁▁▁ |
| numeric | Count_Initial | 0 | 1 | NA | NA | NA | NA | NA | 1.837963 | 2.050129 | 0 | 0.00 | 1.0 | 3.00 | 11 | ▇▂▁▁▁ |
- Creating new columns
Count_germlings$TankCode <- paste(Count_germlings$Treatment,Count_germlings$Tank)
Count_germlings$Well <- paste(Count_germlings$Tank,Count_germlings$Plate,Count_germlings$Number)
#Calculate initial per tank
Tank_Count1 <- Count_germlings%>%
filter(Date == "22/11/2022") %>%
group_by(Site, Treatment, TankCode)%>%
summarise(Sum_count = sum(Count_Initial))
CountingTank <- left_join(Tank_Count1, Count_germlings, by =c('Site'='Site', 'Treatment'='Treatment', 'TankCode'='TankCode'))- Calculate mean Survival by tank per total
Count_sumTank <- CountingTank %>%
group_by(Site, Treatment, TankCode, Day, Sum_count)%>%
summarise(Sum_Tank = sum(Count_Alive))
percentTank <- Count_sumTank %>%
group_by(Site, Treatment, TankCode, Day)%>%
summarise(Percent_survival = ((Sum_Tank*100)/Sum_count)) %>%
filter(Percent_survival != "NaN")
Average_percentTank <- percentTank%>%
group_by(Site, Treatment, Day)%>%
summarise(Percent_mean = mean(Percent_survival),n = n(), sd = sd(Percent_survival), se = sd/sqrt(n)) 2. Germling status
- Check data structure
Length_germlings %>%
skimr::skim() %>%
gt::gt()| skim_type | skim_variable | n_missing | complete_rate | character.min | character.max | character.empty | character.n_unique | character.whitespace | numeric.mean | numeric.sd | numeric.p0 | numeric.p25 | numeric.p50 | numeric.p75 | numeric.p100 | numeric.hist |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| character | Treatment_Code | 0 | 1 | 2 | 2 | 0 | 3 | 0 | NA | NA | NA | NA | NA | NA | NA | NA |
| character | Treatment | 0 | 1 | 4 | 7 | 0 | 3 | 0 | NA | NA | NA | NA | NA | NA | NA | NA |
| character | Plate | 0 | 1 | 1 | 1 | 0 | 2 | 0 | NA | NA | NA | NA | NA | NA | NA | NA |
| character | Site | 0 | 1 | 10 | 10 | 0 | 2 | 0 | NA | NA | NA | NA | NA | NA | NA | NA |
| character | Date | 0 | 1 | 9 | 10 | 0 | 4 | 0 | NA | NA | NA | NA | NA | NA | NA | NA |
| character | Shape | 0 | 1 | 7 | 10 | 0 | 4 | 0 | NA | NA | NA | NA | NA | NA | NA | NA |
| character | Shape2 | 0 | 1 | 7 | 9 | 0 | 3 | 0 | NA | NA | NA | NA | NA | NA | NA | NA |
| character | Status | 0 | 1 | 5 | 9 | 0 | 3 | 0 | NA | NA | NA | NA | NA | NA | NA | NA |
| numeric | Tank | 0 | 1 | NA | NA | NA | NA | NA | 3.3408304 | 1.71437818 | 1.000 | 2.00000 | 3.0000 | 5.00000 | 6.000 | ▇▃▃▃▃ |
| numeric | Number | 0 | 1 | NA | NA | NA | NA | NA | 3.6799308 | 1.70013110 | 1.000 | 2.00000 | 4.0000 | 5.00000 | 6.000 | ▇▅▃▅▅ |
| numeric | Replicate | 0 | 1 | NA | NA | NA | NA | NA | 14.0536332 | 9.09999689 | 1.000 | 6.00000 | 13.0000 | 21.00000 | 37.000 | ▇▆▅▃▂ |
| numeric | Day | 0 | 1 | NA | NA | NA | NA | NA | 19.2335640 | 7.88252887 | 9.000 | 9.00000 | 16.0000 | 23.00000 | 30.000 | ▇▇▁▇▇ |
| numeric | Length | 0 | 1 | NA | NA | NA | NA | NA | 1.1367561 | 0.41711936 | 0.300 | 0.80700 | 1.1165 | 1.41375 | 2.649 | ▅▇▆▂▁ |
| numeric | Length_noriz | 0 | 1 | NA | NA | NA | NA | NA | 0.7819066 | 0.30142459 | 0.204 | 0.54500 | 0.7420 | 0.96375 | 1.737 | ▅▇▆▂▁ |
| numeric | Area | 0 | 1 | NA | NA | NA | NA | NA | 0.1660502 | 0.09626433 | 0.020 | 0.08525 | 0.1470 | 0.22750 | 0.531 | ▇▆▃▁▁ |
- Calculate mean length by tank per total
Length_germlings$TankCode <- paste(Length_germlings$Treatment,Length_germlings$Tank)
#Calculate average per well
Length_tank <- Length_germlings%>%
group_by(Site, Treatment, TankCode, Day)%>%
summarise(Length_pertank = mean(Length_noriz))
Average_Length <- Length_tank%>%
group_by(Site, Treatment, Day)%>%
summarise(MeanLength = mean(Length_pertank),n = n(), sd = sd(Length_pertank), se = sd/sqrt(n))- Calculate mean area by tank per total
#Calculate average per well
Area_tank <- Length_germlings%>%
group_by(Site, Treatment, TankCode, Day)%>%
summarise(Well_meanArea = mean(Area))
Average_Area <- Area_tank%>%
group_by(Site, Treatment, Day)%>%
summarise(MeanArea = mean(Well_meanArea),n = n(), sd = sd(Well_meanArea), se = sd/sqrt(n))3. Germlings shape
- Calculate total numbers for each shape
Length_shapes <- Length_germlings %>%
group_by(Site, Treatment, Day, TankCode, Shape2) %>%
summarise(number = n())
Shape_wide <- Length_shapes %>%
ungroup() %>%
pivot_wider(names_from = Shape2, values_from = number,
values_fill = 0)
Shape_long <- Shape_wide %>%
gather(Shape2, number, Elongated:Divided)
Shape_long <- Shape_long %>%
group_by(Site, Treatment, Day, TankCode) %>%
mutate(percent = number/sum(number), n = sum(number))- Calculate mean per group
Average_Shape <- Shape_long%>%
group_by(Site, Treatment, Day, Shape2)%>%
summarise(MeanShape = mean(percent),n = n(), sd = sd(percent), se = sd/sqrt(n))4. Tanks_metadata
This data correspond to the location of each individual 2 L container in the experiment.
Tanks_metadata %>%
skimr::skim() %>%
gt::gt()| skim_type | skim_variable | n_missing | complete_rate | character.min | character.max | character.empty | character.n_unique | character.whitespace | numeric.mean | numeric.sd | numeric.p0 | numeric.p25 | numeric.p50 | numeric.p75 | numeric.p100 | numeric.hist |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| character | Code | 0 | 1 | 4 | 6 | 0 | 144 | 0 | NA | NA | NA | NA | NA | NA | NA | NA |
| character | Tank_Code | 0 | 1 | 4 | 4 | 0 | 18 | 0 | NA | NA | NA | NA | NA | NA | NA | NA |
| character | Shelve | 0 | 1 | 2 | 4 | 0 | 2 | 0 | NA | NA | NA | NA | NA | NA | NA | NA |
| character | Location | 0 | 1 | 4 | 6 | 0 | 6 | 0 | NA | NA | NA | NA | NA | NA | NA | NA |
| numeric | Rack | 0 | 1 | NA | NA | NA | NA | NA | 2 | 0.8193465 | 1 | 1 | 2 | 3 | 3 | ▇▁▇▁▇ |
5. Juvenile status, survival and length
- Check data structure
Juveniles %>%
skimr::skim() %>%
gt::gt()| skim_type | skim_variable | n_missing | complete_rate | character.min | character.max | character.empty | character.n_unique | character.whitespace | numeric.mean | numeric.sd | numeric.p0 | numeric.p25 | numeric.p50 | numeric.p75 | numeric.p100 | numeric.hist |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| character | Treatment | 0 | 1.0000000 | 4 | 7 | 0 | 3 | 0 | NA | NA | NA | NA | NA | NA | NA | NA |
| character | Site | 0 | 1.0000000 | 8 | 11 | 0 | 4 | 0 | NA | NA | NA | NA | NA | NA | NA | NA |
| character | Distribution | 0 | 1.0000000 | 7 | 9 | 0 | 2 | 0 | NA | NA | NA | NA | NA | NA | NA | NA |
| character | Site_Code | 0 | 1.0000000 | 2 | 2 | 0 | 4 | 0 | NA | NA | NA | NA | NA | NA | NA | NA |
| character | Code | 0 | 1.0000000 | 4 | 6 | 0 | 122 | 0 | NA | NA | NA | NA | NA | NA | NA | NA |
| character | Date | 0 | 1.0000000 | 9 | 10 | 0 | 10 | 0 | NA | NA | NA | NA | NA | NA | NA | NA |
| character | Stage | 0 | 1.0000000 | 5 | 16 | 0 | 9 | 0 | NA | NA | NA | NA | NA | NA | NA | NA |
| character | Measurement | 0 | 1.0000000 | 0 | 7 | 774 | 3 | 0 | NA | NA | NA | NA | NA | NA | NA | NA |
| character | Condition | 0 | 1.0000000 | 0 | 5 | 10 | 9 | 0 | NA | NA | NA | NA | NA | NA | NA | NA |
| character | StipeRoot | 0 | 1.0000000 | 0 | 3 | 85 | 3 | 0 | NA | NA | NA | NA | NA | NA | NA | NA |
| character | Holdfast | 0 | 1.0000000 | 0 | 3 | 86 | 3 | 0 | NA | NA | NA | NA | NA | NA | NA | NA |
| character | Survival | 0 | 1.0000000 | 4 | 5 | 0 | 2 | 0 | NA | NA | NA | NA | NA | NA | NA | NA |
| numeric | Replicate | 0 | 1.0000000 | NA | NA | NA | NA | NA | 6.3543307 | 3.5023484 | 1.00 | 3.000 | 6.000 | 9.00 | 14.000 | ▇▇▅▇▂ |
| numeric | Year | 0 | 1.0000000 | NA | NA | NA | NA | NA | 2023.0000000 | 0.0000000 | 2023.00 | 2023.000 | 2023.000 | 2023.00 | 2023.000 | ▁▁▇▁▁ |
| numeric | Day_equivalent | 0 | 1.0000000 | NA | NA | NA | NA | NA | 14.3110236 | 7.0503968 | 0.00 | 12.000 | 16.000 | 20.00 | 25.000 | ▂▂▅▇▃ |
| numeric | Length_cm | 683 | 0.3277559 | NA | NA | NA | NA | NA | 12.0765766 | 4.7646974 | 1.00 | 9.000 | 11.700 | 14.80 | 32.000 | ▂▇▃▁▁ |
| numeric | Initial_thickness_mm | 50 | 0.9507874 | NA | NA | NA | NA | NA | 0.6102381 | 0.1198775 | 0.37 | 0.520 | 0.590 | 0.68 | 0.990 | ▂▇▅▂▁ |
| numeric | Fv_Fm | 583 | 0.4261811 | NA | NA | NA | NA | NA | 0.3167067 | 0.1534231 | 0.00 | 0.238 | 0.312 | 0.40 | 0.875 | ▂▇▅▁▁ |
| numeric | Melting | 85 | 0.9163386 | NA | NA | NA | NA | NA | 0.9849624 | 1.2842258 | 0.00 | 0.000 | 1.000 | 1.00 | 4.000 | ▇▅▁▂▂ |
| numeric | Bleaching | 86 | 0.9153543 | NA | NA | NA | NA | NA | 0.2903226 | 0.5099789 | 0.00 | 0.000 | 0.000 | 1.00 | 3.000 | ▇▂▁▁▁ |
| numeric | Fouling | 86 | 0.9153543 | NA | NA | NA | NA | NA | 0.3806452 | 0.5560622 | 0.00 | 0.000 | 0.000 | 1.00 | 2.000 | ▇▁▃▁▁ |
| numeric | Black | 86 | 0.9153543 | NA | NA | NA | NA | NA | 0.4204301 | 0.7399210 | 0.00 | 0.000 | 0.000 | 1.00 | 4.000 | ▇▂▁▁▁ |
| numeric | Green | 87 | 0.9143701 | NA | NA | NA | NA | NA | 0.1959096 | 0.5554825 | 0.00 | 0.000 | 0.000 | 0.00 | 4.000 | ▇▁▁▁▁ |
| numeric | N_blades | 698 | 0.3129921 | NA | NA | NA | NA | NA | 4.3899371 | 1.7143372 | 0.00 | 3.000 | 4.000 | 6.00 | 11.000 | ▂▇▇▁▁ |
| numeric | Survival_code | 0 | 1.0000000 | NA | NA | NA | NA | NA | 0.8444882 | 0.3625704 | 0.00 | 1.000 | 1.000 | 1.00 | 1.000 | ▂▁▁▁▇ |
- Join juvenile measurements and metadata
Juveniles_meta <- Juveniles %>%
left_join(Tanks_metadata, by = c("Code")) - Calculate means for melting, bleaching, fouling, black + green
Juv_means <- Juveniles_meta %>%
group_by(Distribution, Site, Treatment, Day_equivalent) %>%
summarise(Mean_melting = mean(Melting, na.rm=T),n = n(), sd_melting = sd(Melting, na.rm=T), se_melting = sd_melting/sqrt(n),
Mean_bleaching = mean(Bleaching, na.rm=T), sd_bleaching = sd(Bleaching, na.rm=T), se_bleaching = sd_bleaching/sqrt(n),
Mean_fouling = mean(Fouling, na.rm=T), sd_fouling = sd(Fouling, na.rm=T), se_fouling = sd_fouling/sqrt(n),
Mean_necrosis = mean(Black, na.rm=T), sd_necrosis = sd(Black, na.rm=T), se_necrosis = sd_necrosis/sqrt(n),
Mean_green = mean(Green, na.rm=T), sd_green = sd(Green, na.rm=T), se_green = sd_green/sqrt(n),
Mean_survival = mean(Survival_code, na.rm=T), sd_surv = sd(Survival_code, na.rm=T), se_surv = sd_surv/sqrt(n)) - Calculate means for Fv/Fm
Photo_means <- Juveniles_meta %>%
filter(Day_equivalent == 8| Day_equivalent == 12| Day_equivalent == 18| Day_equivalent == 22) %>%
group_by(Distribution, Site, Treatment, Day_equivalent) %>%
summarise(Mean_FvFm = mean(Fv_Fm, na.rm=T),n = n(), sd_FvFm = sd(Fv_Fm, na.rm=T), se_FvFm = sd_FvFm/sqrt(n),) - Calculate survival (%)
#Count initial numbers per treatment and day
Juv_count <- Juveniles_meta%>%
filter(Day_equivalent == "0") %>%
group_by(Site, Treatment)%>%
summarise(Initial_count = n())
Counting <- left_join(Juv_count, Juveniles_meta, by =c('Site'='Site', 'Treatment'='Treatment'))
Juv_count2 <- Counting %>%
group_by(Site, Distribution, Tank_Code, Treatment, Day_equivalent, Initial_count)%>%
summarise(Alive_juv = sum(Survival == "Alive", na.rm=TRUE))
Percent_Survival <- Juv_count2 %>%
group_by(Site, Distribution, Tank_Code, Treatment, Day_equivalent)%>%
summarise(Percent_survival = ((Alive_juv*100)/Initial_count))
##Calculate summary data
Survival_means <- Percent_Survival %>%
filter(Day_equivalent == 8| Day_equivalent == 12| Day_equivalent == 18| Day_equivalent == 22) %>%
group_by(Distribution, Treatment, Day_equivalent) %>%
summarise(Mean_survival = mean(Percent_survival, na.rm=T),n = n(), sd = sd(Percent_survival, na.rm=T), se = sd/sqrt(n)) - Calculate length means
Juveniles_length <- Juveniles_meta %>%
filter(Measurement == "Initial" | Measurement == "Final") %>%
filter(Survival != "Dead") %>%
group_by(Code, Treatment, Distribution, Site, Tank_Code, Location) %>%
summarise(diff = (((Length_cm[match("Final", Measurement)]/Length_cm[match("Initial", Measurement)])^(1/17))-1)*100)
Length_means <- Juveniles_length %>%
group_by(Site, Treatment) %>%
summarise(Mean_length = mean(diff, na.rm=T),n = n(), sd = sd(diff, na.rm=T), se = sd/sqrt(n)) 6. Juvenile CN results
- Check data structure
CN_results %>%
skimr::skim() %>%
gt::gt()| skim_type | skim_variable | n_missing | complete_rate | character.min | character.max | character.empty | character.n_unique | character.whitespace | numeric.mean | numeric.sd | numeric.p0 | numeric.p25 | numeric.p50 | numeric.p75 | numeric.p100 | numeric.hist |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| character | Treatment | 0 | 1 | 4 | 7 | 0 | 3 | 0 | NA | NA | NA | NA | NA | NA | NA | NA |
| character | Site | 0 | 1 | 8 | 11 | 0 | 4 | 0 | NA | NA | NA | NA | NA | NA | NA | NA |
| character | Code | 0 | 1 | 4 | 6 | 0 | 71 | 0 | NA | NA | NA | NA | NA | NA | NA | NA |
| character | Method | 0 | 1 | 14 | 14 | 0 | 1 | 0 | NA | NA | NA | NA | NA | NA | NA | NA |
| character | Name | 0 | 1 | 1 | 6 | 0 | 75 | 0 | NA | NA | NA | NA | NA | NA | NA | NA |
| character | Sample.Code | 0 | 1 | 4 | 10 | 0 | 74 | 0 | NA | NA | NA | NA | NA | NA | NA | NA |
| numeric | Number | 0 | 1 | NA | NA | NA | NA | NA | 38.053333 | 2.188353e+01 | 1.0000 | 19.5000 | 38.0000 | 56.5000 | 76.0000 | ▇▇▇▇▇ |
| numeric | Weight_mg | 0 | 1 | NA | NA | NA | NA | NA | 2.064643 | 9.140592e-02 | 1.8420 | 2.0075 | 2.0600 | 2.1260 | 2.3640 | ▂▇▇▂▁ |
| numeric | N.Area | 0 | 1 | NA | NA | NA | NA | NA | 587.346667 | 1.178839e+02 | 382.0000 | 502.0000 | 567.0000 | 665.5000 | 958.0000 | ▅▇▆▂▁ |
| numeric | C.Area | 0 | 1 | NA | NA | NA | NA | NA | 15173.826667 | 1.387957e+03 | 11758.0000 | 14167.0000 | 15249.0000 | 16175.5000 | 17636.0000 | ▂▅▇▇▅ |
| numeric | N | 0 | 1 | NA | NA | NA | NA | NA | 0.880000 | 1.663113e-01 | 0.5700 | 0.7700 | 0.8500 | 0.9900 | 1.4300 | ▅▇▇▂▁ |
| numeric | C | 0 | 1 | NA | NA | NA | NA | NA | 31.915733 | 2.433816e+00 | 25.8600 | 30.4800 | 32.0000 | 33.7600 | 36.8600 | ▂▃▇▆▂ |
| numeric | CN_ratio | 0 | 1 | NA | NA | NA | NA | NA | 37.657040 | 8.261379e+00 | 22.0217 | 31.6244 | 36.4972 | 41.1365 | 58.4604 | ▃▇▆▃▂ |
| numeric | N.Factor | 0 | 1 | NA | NA | NA | NA | NA | 1.071600 | 0.000000e+00 | 1.0716 | 1.0716 | 1.0716 | 1.0716 | 1.0716 | ▁▁▇▁▁ |
| numeric | C.Factor | 0 | 1 | NA | NA | NA | NA | NA | 1.085400 | 0.000000e+00 | 1.0854 | 1.0854 | 1.0854 | 1.0854 | 1.0854 | ▁▁▇▁▁ |
| numeric | N.Blank | 0 | 1 | NA | NA | NA | NA | NA | 8.000000 | 0.000000e+00 | 8.0000 | 8.0000 | 8.0000 | 8.0000 | 8.0000 | ▁▁▇▁▁ |
| numeric | C.Blank | 0 | 1 | NA | NA | NA | NA | NA | 93.000000 | 0.000000e+00 | 93.0000 | 93.0000 | 93.0000 | 93.0000 | 93.0000 | ▁▁▇▁▁ |
- Join CN data with individual characteristics
Juv_1 <- filter(Juveniles, Measurement == "Final")
CN_results1 <- CN_results %>%
group_by(Treatment, Site, Code)%>%
summarise(CN_ratio = mean(CN_ratio), C = mean(C), N = mean(N))
CN_results2 <- inner_join(CN_results1, Juv_1, by=c("Code"))- Maintain only alive individuals from MHW1 and ambient
CN_results3 <- filter(CN_results2, Survival !="DEAD")
CN_results4 <- filter(CN_results3, Treatment.y !="MHW2")- Calculate average values
CNresults_mean <- CN_results4 %>%
group_by(Treatment.y, Site.x) %>%
summarise(Mean_CNratio = mean(CN_ratio, na.rm=T),n = n(), sd = sd(CN_ratio, na.rm=T), se = sd/sqrt(n), Mean_C = mean(C, na.rm=T),n = n(), sd_C = sd(C, na.rm=T), se_C = sd_C/sqrt(n), Mean_N = mean(N, na.rm=T),n = n(), sd_N = sd(N, na.rm=T), se_N = sd_N/sqrt(n))7. Juvenile pigment results
- Join pigment data and check data structure
Pigments2 <- inner_join(Pigments, Acetone, by=c("Sample"))
Pigments2 <- inner_join(Pigments2, DMSO, by=c("Sample"))
Pigments2 %>%
skimr::skim() %>%
gt::gt()| skim_type | skim_variable | n_missing | complete_rate | character.min | character.max | character.empty | character.n_unique | character.whitespace | numeric.mean | numeric.sd | numeric.p0 | numeric.p25 | numeric.p50 | numeric.p75 | numeric.p100 | numeric.hist |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| character | Treatment | 0 | 1 | 4 | 7 | 0 | 4 | 0 | NA | NA | NA | NA | NA | NA | NA | NA |
| character | Site_Code | 0 | 1 | 8 | 11 | 0 | 4 | 0 | NA | NA | NA | NA | NA | NA | NA | NA |
| character | Code | 0 | 1 | 4 | 6 | 0 | 92 | 0 | NA | NA | NA | NA | NA | NA | NA | NA |
| numeric | Sample | 0 | 1 | NA | NA | NA | NA | NA | 57.333333333 | 3.387801e+01 | 1.000000000 | 2.875000e+01 | 55.500000000 | 87.250000000 | 1.150000e+02 | ▇▇▇▇▇ |
| numeric | Weight | 0 | 1 | NA | NA | NA | NA | NA | 0.139234259 | 1.340195e-02 | 0.091100000 | 1.324500e-01 | 0.144100000 | 0.149525000 | 1.597000e-01 | ▁▁▃▅▇ |
| numeric | X470.0nm | 0 | 1 | NA | NA | NA | NA | NA | 0.498842593 | 2.712479e-01 | 0.134800000 | 3.008250e-01 | 0.442500000 | 0.678025000 | 1.560500e+00 | ▇▅▃▁▁ |
| numeric | X581.0nm | 0 | 1 | NA | NA | NA | NA | NA | 0.135269444 | 5.958827e-02 | 0.061800000 | 9.312500e-02 | 0.124200000 | 0.157725000 | 4.048000e-01 | ▇▅▂▁▁ |
| numeric | X631.0nm.x | 0 | 1 | NA | NA | NA | NA | NA | 0.154786111 | 7.192443e-02 | 0.065700000 | 1.007000e-01 | 0.143650000 | 0.182150000 | 4.722000e-01 | ▇▆▂▁▁ |
| numeric | X664.0nm | 0 | 1 | NA | NA | NA | NA | NA | 0.572501852 | 3.268953e-01 | 0.145600000 | 3.246000e-01 | 0.514850000 | 0.739225000 | 1.895400e+00 | ▇▆▂▁▁ |
| numeric | X750.0nm.x | 0 | 1 | NA | NA | NA | NA | NA | 0.051438889 | 8.732717e-03 | 0.035300000 | 4.627500e-02 | 0.049200000 | 0.053650000 | 9.890000e-02 | ▆▇▂▁▁ |
| numeric | Chla_acetone | 0 | 1 | NA | NA | NA | NA | NA | 0.007778558 | 4.441512e-03 | 0.001978261 | 4.410327e-03 | 0.006995245 | 0.010043818 | 2.575272e-02 | ▇▆▂▁▁ |
| numeric | Chlc_acetone | 0 | 1 | NA | NA | NA | NA | NA | 0.001902010 | 6.575584e-04 | 0.001131672 | 1.461254e-03 | 0.001681993 | 0.002124558 | 4.957878e-03 | ▇▃▁▁▁ |
| numeric | Fucoxanthin_acetone | 0 | 1 | NA | NA | NA | NA | NA | 0.002386659 | 1.652937e-03 | -0.000028600 | 1.135309e-03 | 0.002013204 | 0.003409626 | 8.693583e-03 | ▇▆▃▁▁ |
| numeric | X480.0nm | 0 | 1 | NA | NA | NA | NA | NA | 0.397030556 | 3.802401e-01 | 0.090200000 | 1.566500e-01 | 0.223700000 | 0.488525000 | 2.119700e+00 | ▇▁▁▁▁ |
| numeric | X582.0nm | 0 | 1 | NA | NA | NA | NA | NA | 0.131809259 | 1.033201e-01 | 0.052800000 | 7.170000e-02 | 0.085850000 | 0.136350000 | 5.340000e-01 | ▇▁▁▁▁ |
| numeric | X631.0nm.y | 0 | 1 | NA | NA | NA | NA | NA | 0.152635185 | 1.342896e-01 | 0.053900000 | 7.430000e-02 | 0.091300000 | 0.158200000 | 6.961000e-01 | ▇▁▁▁▁ |
| numeric | X665.0nm | 0 | 1 | NA | NA | NA | NA | NA | 0.456497222 | 5.175381e-01 | 0.078900000 | 1.511250e-01 | 0.219100000 | 0.463775000 | 2.448400e+00 | ▇▁▁▁▁ |
| numeric | X750.0nm.y | 0 | 1 | NA | NA | NA | NA | NA | 0.055616667 | 2.011757e-02 | 0.040000000 | 4.460000e-02 | 0.048200000 | 0.055000000 | 1.298000e-01 | ▇▁▁▁▁ |
| numeric | Chla_DMSO | 0 | 1 | NA | NA | NA | NA | NA | 0.006270566 | 7.109040e-03 | 0.001083791 | 2.075893e-03 | 0.003009615 | 0.006370536 | 3.363187e-02 | ▇▁▁▁▁ |
| numeric | Chlc_DMSO | 0 | 1 | NA | NA | NA | NA | NA | 0.002408815 | 1.389306e-03 | 0.001249799 | 1.581211e-03 | 0.001807062 | 0.002523921 | 8.303749e-03 | ▇▁▁▁▁ |
| numeric | Fucoxanthin_DMSO | 0 | 1 | NA | NA | NA | NA | NA | 0.002055245 | 2.309433e-03 | 0.000136831 | 5.807765e-04 | 0.001022804 | 0.002522441 | 1.283262e-02 | ▇▁▁▁▁ |
- Calculate means
Pigments3 <- Pigments2 %>%
group_by(Treatment, Site_Code, Code)%>%
summarise(Chla_acetone = mean(Chla_acetone), Chlc_acetone = mean(Chlc_acetone), Fucoxanthin_acetone = mean(Fucoxanthin_acetone), Weight = mean(Weight), Chla_DMSO = mean(Chla_DMSO), Chlc_DMSO= mean(Chlc_DMSO), Fucoxanthin_DMSO = mean(Fucoxanthin_DMSO))- Calculate total pigment content
Pigments3$chla <- Pigments3$Chla_acetone +Pigments3$Chla_DMSO
Pigments3$chlc <- Pigments3$Chlc_acetone +Pigments3$Chlc_DMSO
Pigments3$Fucoxanthin <- Pigments3$Fucoxanthin_acetone +Pigments3$Fucoxanthin_DMSO- Mantain data from alive individuals from MHW1 and ambient
Pigments4 <- inner_join(Pigments3, Juv_1, by=c("Code"))
Pigments5 <- filter(Pigments4, Survival != "DEAD")
Pigments6 <- filter(Pigments5, Treatment.y != "MHW2")- Calculate average values
Pigments_mean <- Pigments6 %>%
group_by(Treatment.y, Site_Code.x) %>%
summarise(Mean_chla = mean(chla, na.rm=T),n = n(), sd_chla = sd(chla, na.rm=T), se_chla = sd_chla/sqrt(n), Mean_Fucoxanthin = mean(Fucoxanthin, na.rm=T),n = n(), sd_Fucoxanthin = sd(Fucoxanthin, na.rm=T), se_Fucoxanthin = sd_Fucoxanthin/sqrt(n), Mean_chlc = mean(chlc, na.rm=T),n = n(), sd_chlc = sd(chlc, na.rm=T), se_chlc = sd_chlc/sqrt(n))Section 2: Germling experiments
Data from germling experimnets with individuals from Black Head and Palm Beach was visualised and analised to evaluate differences among populations and MHW treatments
Set general theme for plots
Theme1 <- theme_bw() +
theme(panel.grid.major = element_line(colour = "white"),
panel.grid.minor = element_line(colour = "white"), axis.text.y = element_text(size = 12), axis.title.y = element_text(size = 12), axis.text.x = element_text(size = 12), axis.title.x = element_text(size = 12), legend.title = element_text(size=12), #change legend title font size
legend.text = element_text(size=12))2a) Germling survival
1. Visualise differences
Average_percentTank2 <- Average_percentTank
Average_percentTank2$Site <- factor(
Average_percentTank2$Site,
levels = c("Black_Head", "Palm_Beach"),
labels = c(
"Black Head", "Palm Beach"
))Create plot to visualise differences
Average_percentTank2$Day <- as.numeric(Average_percentTank2$Day)
Count_plot16 <- Average_percentTank2 %>%
filter(Site == "Palm Beach" | Site == "Black Head") %>%
filter(Day == 9|Day == 16|Day == 23| Day == 30) %>%
ggplot(aes(x= Day, y= Percent_mean, colour = Site)) +
geom_errorbar(aes(ymin=Percent_mean-se, ymax=Percent_mean+se),colour = "grey", width=.2) +
geom_point() +
geom_line()+
Theme1 +
labs( y="Survival (%)", x="Time (days)") +
scale_colour_manual(values = c("orange", "darkolivegreen3")) +
theme(strip.text.x = element_text(size = 12))+
scale_y_continuous(limits = c(0, 105), breaks = seq(0, 100, by = 25))+
scale_x_continuous(limits = c(8, 31), breaks = c(9, 16, 23,30))+
facet_wrap(~ factor(Treatment, c("Ambient","MHW1","MHW2")), ncol = 3)
Count_plot16
ggsave("Survival_germlings.jpeg", width = 200, height = 56, units = "mm")2. Statistical analyses
Filter data per day
Percent_Tank <- filter (percentTank, Site == "Palm_Beach" | Site == "Black_Head")
Percent_Tank$Day <- as.factor(Percent_Tank$Day)
Percent_Tank$Proportion_Survival <- (Percent_Tank$Percent_survival)/100
Percent_TankALL <- filter (Percent_Tank, Day == "9" |Day == "16"|Day == "23" |Day == "30")- Run LMM for germling survival
my_lm1 <- lmer(Proportion_Survival ~ Site*Treatment*Day + (1|TankCode), data=Percent_TankALL)
summary(my_lm1)## Linear mixed model fit by REML ['lmerMod']
## Formula: Proportion_Survival ~ Site * Treatment * Day + (1 | TankCode)
## Data: Percent_TankALL
##
## REML criterion at convergence: -48.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.1660 -0.3941 0.0638 0.3180 3.4963
##
## Random effects:
## Groups Name Variance Std.Dev.
## TankCode (Intercept) 0.009781 0.0989
## Residual 0.022633 0.1504
## Number of obs: 144, groups: TankCode, 18
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.000e+00 7.350e-02 13.605
## SitePalm_Beach -1.925e-16 8.686e-02 0.000
## TreatmentMHW1 -3.617e-15 1.039e-01 0.000
## TreatmentMHW2 -5.550e-15 1.039e-01 0.000
## Day16 -2.381e-02 8.686e-02 -0.274
## Day23 -5.159e-02 8.686e-02 -0.594
## Day30 -1.488e-01 8.686e-02 -1.713
## SitePalm_Beach:TreatmentMHW1 -3.371e-15 1.228e-01 0.000
## SitePalm_Beach:TreatmentMHW2 -6.119e-16 1.228e-01 0.000
## SitePalm_Beach:Day16 2.381e-02 1.228e-01 0.194
## SitePalm_Beach:Day23 2.937e-02 1.228e-01 0.239
## SitePalm_Beach:Day30 9.464e-02 1.228e-01 0.770
## TreatmentMHW1:Day16 -4.286e-02 1.228e-01 -0.349
## TreatmentMHW2:Day16 -3.623e-01 1.228e-01 -2.949
## TreatmentMHW1:Day23 -1.508e-02 1.228e-01 -0.123
## TreatmentMHW2:Day23 -6.123e-01 1.228e-01 -4.985
## TreatmentMHW1:Day30 4.881e-02 1.228e-01 0.397
## TreatmentMHW2:Day30 -5.151e-01 1.228e-01 -4.193
## SitePalm_Beach:TreatmentMHW1:Day16 3.452e-02 1.737e-01 0.199
## SitePalm_Beach:TreatmentMHW2:Day16 -1.345e-01 1.737e-01 -0.774
## SitePalm_Beach:TreatmentMHW1:Day23 -4.980e-02 1.737e-01 -0.287
## SitePalm_Beach:TreatmentMHW2:Day23 -1.266e-01 1.737e-01 -0.729
## SitePalm_Beach:TreatmentMHW1:Day30 -8.174e-02 1.737e-01 -0.471
## SitePalm_Beach:TreatmentMHW2:Day30 -2.268e-01 1.737e-01 -1.305##
## Correlation matrix not shown by default, as p = 24 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need itAnova(my_lm1)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Proportion_Survival
## Chisq Df Pr(>Chisq)
## Site 0.2210 1 0.6383
## Treatment 57.6670 2 3.004e-13 ***
## Day 92.3452 3 < 2.2e-16 ***
## Site:Treatment 4.4208 2 0.1097
## Site:Day 0.1862 3 0.9798
## Treatment:Day 98.2682 6 < 2.2e-16 ***
## Site:Treatment:Day 2.2690 6 0.8934
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1emmeans(my_lm1, specs = pairwise ~ Treatment|Day, type = "response", adjust = "none")## NOTE: Results may be misleading due to involvement in interactions## $emmeans
## Day = 9:
## Treatment emmean SE df lower.CL upper.CL
## Ambient 1.000 0.0593 39.4 0.880 1.120
## MHW1 1.000 0.0593 39.4 0.880 1.120
## MHW2 1.000 0.0593 39.4 0.880 1.120
##
## Day = 16:
## Treatment emmean SE df lower.CL upper.CL
## Ambient 0.988 0.0593 39.4 0.868 1.108
## MHW1 0.963 0.0593 39.4 0.843 1.082
## MHW2 0.559 0.0593 39.4 0.439 0.678
##
## Day = 23:
## Treatment emmean SE df lower.CL upper.CL
## Ambient 0.963 0.0593 39.4 0.843 1.083
## MHW1 0.923 0.0593 39.4 0.803 1.043
## MHW2 0.287 0.0593 39.4 0.168 0.407
##
## Day = 30:
## Treatment emmean SE df lower.CL upper.CL
## Ambient 0.899 0.0593 39.4 0.779 1.018
## MHW1 0.906 0.0593 39.4 0.787 1.026
## MHW2 0.270 0.0593 39.4 0.150 0.390
##
## Results are averaged over the levels of: Site
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## Day = 9:
## contrast estimate SE df t.ratio p.value
## Ambient - MHW1 0.00000 0.0839 39.4 0.000 1.0000
## Ambient - MHW2 0.00000 0.0839 39.4 0.000 1.0000
## MHW1 - MHW2 0.00000 0.0839 39.4 0.000 1.0000
##
## Day = 16:
## contrast estimate SE df t.ratio p.value
## Ambient - MHW1 0.02560 0.0839 39.4 0.305 0.7618
## Ambient - MHW2 0.42956 0.0839 39.4 5.122 <.0001
## MHW1 - MHW2 0.40396 0.0839 39.4 4.817 <.0001
##
## Day = 23:
## contrast estimate SE df t.ratio p.value
## Ambient - MHW1 0.03998 0.0839 39.4 0.477 0.6362
## Ambient - MHW2 0.67560 0.0839 39.4 8.056 <.0001
## MHW1 - MHW2 0.63562 0.0839 39.4 7.579 <.0001
##
## Day = 30:
## contrast estimate SE df t.ratio p.value
## Ambient - MHW1 -0.00794 0.0839 39.4 -0.095 0.9251
## Ambient - MHW2 0.62847 0.0839 39.4 7.494 <.0001
## MHW1 - MHW2 0.63641 0.0839 39.4 7.589 <.0001
##
## Results are averaged over the levels of: Site
## Degrees-of-freedom method: kenward-rogeremmeans(my_lm1, specs = pairwise ~ Day|Treatment, type = "response", adjust = "none")## NOTE: Results may be misleading due to involvement in interactions## $emmeans
## Treatment = Ambient:
## Day emmean SE df lower.CL upper.CL
## 9 1.000 0.0593 39.4 0.880 1.120
## 16 0.988 0.0593 39.4 0.868 1.108
## 23 0.963 0.0593 39.4 0.843 1.083
## 30 0.899 0.0593 39.4 0.779 1.018
##
## Treatment = MHW1:
## Day emmean SE df lower.CL upper.CL
## 9 1.000 0.0593 39.4 0.880 1.120
## 16 0.963 0.0593 39.4 0.843 1.082
## 23 0.923 0.0593 39.4 0.803 1.043
## 30 0.906 0.0593 39.4 0.787 1.026
##
## Treatment = MHW2:
## Day emmean SE df lower.CL upper.CL
## 9 1.000 0.0593 39.4 0.880 1.120
## 16 0.559 0.0593 39.4 0.439 0.678
## 23 0.287 0.0593 39.4 0.168 0.407
## 30 0.270 0.0593 39.4 0.150 0.390
##
## Results are averaged over the levels of: Site
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## Treatment = Ambient:
## contrast estimate SE df t.ratio p.value
## Day9 - Day16 0.0119 0.0614 105 0.194 0.8467
## Day9 - Day23 0.0369 0.0614 105 0.601 0.5492
## Day9 - Day30 0.1015 0.0614 105 1.652 0.1014
## Day16 - Day23 0.0250 0.0614 105 0.407 0.6848
## Day16 - Day30 0.0896 0.0614 105 1.459 0.1477
## Day23 - Day30 0.0646 0.0614 105 1.052 0.2954
##
## Treatment = MHW1:
## contrast estimate SE df t.ratio p.value
## Day9 - Day16 0.0375 0.0614 105 0.611 0.5428
## Day9 - Day23 0.0769 0.0614 105 1.252 0.2134
## Day9 - Day30 0.0935 0.0614 105 1.523 0.1307
## Day16 - Day23 0.0394 0.0614 105 0.641 0.5228
## Day16 - Day30 0.0560 0.0614 105 0.913 0.3635
## Day23 - Day30 0.0167 0.0614 105 0.271 0.7866
##
## Treatment = MHW2:
## contrast estimate SE df t.ratio p.value
## Day9 - Day16 0.4415 0.0614 105 7.188 <.0001
## Day9 - Day23 0.7125 0.0614 105 11.601 <.0001
## Day9 - Day30 0.7300 0.0614 105 11.885 <.0001
## Day16 - Day23 0.2710 0.0614 105 4.413 <.0001
## Day16 - Day30 0.2885 0.0614 105 4.697 <.0001
## Day23 - Day30 0.0175 0.0614 105 0.284 0.7768
##
## Results are averaged over the levels of: Site
## Degrees-of-freedom method: kenward-rogerModel check
res <- resid(my_lm1)
#produce residual vs. fitted plot
plot(fitted(my_lm1), res)
#add a horizontal line at 0
abline(0,0)
#create Q-Q plot for residuals
qqnorm(res)
#add a straight diagonal line to the plot
qqline(res) 
2b) Germlings length
1. Visualise differences
Average_Length2 <- Average_Length
Average_Length2$Site <- factor(
Average_Length2$Site,
levels = c("Black_Head", "Palm_Beach"),
labels = c(
"Black Head", "Palm Beach"
))Average_Length2$Day <- as.numeric(Average_Length2$Day)
Length_plot16 <- Average_Length2 %>%
filter(Site == "Palm Beach" | Site == "Black Head") %>%
filter(Day == 9|Day == 16|Day == 23| Day == 30) %>%
ggplot(aes(x= Day, y= MeanLength, colour = Site)) +
geom_errorbar(aes(ymin=MeanLength-se, ymax=MeanLength+se),colour = "grey", width=.2) +
geom_point() +
geom_line()+
Theme1 +
labs( y="Length (mm)", x="Time (days)") +
scale_colour_manual(values = c("orange", "darkolivegreen3")) +
theme(strip.text.x = element_text(size = 12))+
scale_y_continuous(limits = c(0, 1.4), breaks = seq(0, 1.2, by = 0.4))+
scale_x_continuous(limits = c(8, 31), breaks = c(9, 16, 23,30))+
facet_wrap(~ factor(Treatment, c("Ambient","MHW1","MHW2")), ncol = 3)
Length_plot16
ggsave("Length_germlings.jpeg", width = 200, height = 56, units = "mm")2. Statistical analyses
- Filter days to analyse and transform to factors
Length_tank1 <- filter (Length_tank, Site == "Palm_Beach" | Site == "Black_Head")
Length_tank1$Day <- as.factor(Length_tank1$Day)
Length_wellALL <- filter (Length_tank1, Day == "9" |Day == "16"|Day == "23" |Day == "30")- Run linear model with site, MHW treatment and day as crossed fixed factors and water bath (tank) as random effect.
my_lm2 <- lmer(Length_pertank ~ Site *Treatment*Day + (1|TankCode), data=Length_wellALL)
summary(my_lm2)## Linear mixed model fit by REML ['lmerMod']
## Formula: Length_pertank ~ Site * Treatment * Day + (1 | TankCode)
## Data: Length_wellALL
##
## REML criterion at convergence: -125.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.2558 -0.5199 -0.0870 0.4616 2.5131
##
## Random effects:
## Groups Name Variance Std.Dev.
## TankCode (Intercept) 0.008366 0.09147
## Residual 0.009279 0.09633
## Number of obs: 130, groups: TankCode, 18
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.508597 0.054230 9.379
## SitePalm_Beach 0.004407 0.055614 0.079
## TreatmentMHW1 -0.005136 0.076692 -0.067
## TreatmentMHW2 0.035769 0.076692 0.466
## Day16 0.225199 0.055614 4.049
## Day23 0.324519 0.055614 5.835
## Day30 0.438128 0.055614 7.878
## SitePalm_Beach:TreatmentMHW1 -0.048668 0.078650 -0.619
## SitePalm_Beach:TreatmentMHW2 -0.078440 0.078650 -0.997
## SitePalm_Beach:Day16 -0.045828 0.078650 -0.583
## SitePalm_Beach:Day23 0.006854 0.078650 0.087
## SitePalm_Beach:Day30 0.029225 0.078650 0.372
## TreatmentMHW1:Day16 0.124354 0.078650 1.581
## TreatmentMHW2:Day16 -0.117307 0.084381 -1.390
## TreatmentMHW1:Day23 0.227569 0.078650 2.893
## TreatmentMHW2:Day23 -0.048874 0.084435 -0.579
## TreatmentMHW1:Day30 0.301803 0.078650 3.837
## TreatmentMHW2:Day30 0.060892 0.084435 0.721
## SitePalm_Beach:TreatmentMHW1:Day16 -0.022881 0.111227 -0.206
## SitePalm_Beach:TreatmentMHW2:Day16 0.118281 0.118740 0.996
## SitePalm_Beach:TreatmentMHW1:Day23 -0.131321 0.111227 -1.181
## SitePalm_Beach:TreatmentMHW2:Day23 -0.039591 0.122805 -0.322
## SitePalm_Beach:TreatmentMHW1:Day30 -0.174905 0.111227 -1.573
## SitePalm_Beach:TreatmentMHW2:Day30 -0.286736 0.122805 -2.335##
## Correlation matrix not shown by default, as p = 24 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need itAnova(my_lm2)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Length_pertank
## Chisq Df Pr(>Chisq)
## Site 19.3461 1 1.090e-05 ***
## Treatment 5.2572 2 0.072180 .
## Day 502.4284 3 < 2.2e-16 ***
## Site:Treatment 13.1551 2 0.001391 **
## Site:Day 4.8877 3 0.180210
## Treatment:Day 29.2349 6 5.492e-05 ***
## Site:Treatment:Day 11.7962 6 0.066672 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- Post-hoc test
emmeans(my_lm2, specs = pairwise ~ Treatment|Site, type = "response", adjust = "none")## NOTE: Results may be misleading due to involvement in interactions## $emmeans
## Site = Black_Head:
## Treatment emmean SE df lower.CL upper.CL
## Ambient 0.756 0.0422 17.8 0.667 0.844
## MHW1 0.914 0.0422 17.8 0.825 1.003
## MHW2 0.765 0.0448 21.6 0.672 0.858
##
## Site = Palm_Beach:
## Treatment emmean SE df lower.CL upper.CL
## Ambient 0.758 0.0422 17.8 0.669 0.846
## MHW1 0.785 0.0422 17.8 0.696 0.874
## MHW2 0.637 0.0465 24.6 0.541 0.732
##
## Results are averaged over the levels of: Day
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## Site = Black_Head:
## contrast estimate SE df t.ratio p.value
## Ambient - MHW1 -0.15830 0.0597 17.8 -2.652 0.0163
## Ambient - MHW2 -0.00945 0.0615 19.7 -0.154 0.8796
## MHW1 - MHW2 0.14885 0.0615 19.7 2.419 0.0254
##
## Site = Palm_Beach:
## contrast estimate SE df t.ratio p.value
## Ambient - MHW1 -0.02735 0.0597 17.8 -0.458 0.6523
## Ambient - MHW2 0.12100 0.0628 21.2 1.926 0.0676
## MHW1 - MHW2 0.14835 0.0628 21.2 2.362 0.0279
##
## Results are averaged over the levels of: Day
## Degrees-of-freedom method: kenward-rogeremmeans(my_lm2, specs = pairwise ~ Site|Treatment, type = "response", adjust = "none")## NOTE: Results may be misleading due to involvement in interactions## $emmeans
## Treatment = Ambient:
## Site emmean SE df lower.CL upper.CL
## Black_Head 0.756 0.0422 17.8 0.667 0.844
## Palm_Beach 0.758 0.0422 17.8 0.669 0.846
##
## Treatment = MHW1:
## Site emmean SE df lower.CL upper.CL
## Black_Head 0.914 0.0422 17.8 0.825 1.003
## Palm_Beach 0.785 0.0422 17.8 0.696 0.874
##
## Treatment = MHW2:
## Site emmean SE df lower.CL upper.CL
## Black_Head 0.765 0.0448 21.6 0.672 0.858
## Palm_Beach 0.637 0.0465 24.6 0.541 0.732
##
## Results are averaged over the levels of: Day
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## Treatment = Ambient:
## contrast estimate SE df t.ratio p.value
## Black_Head - Palm_Beach -0.00197 0.0278 91.1 -0.071 0.9437
##
## Treatment = MHW1:
## contrast estimate SE df t.ratio p.value
## Black_Head - Palm_Beach 0.12898 0.0278 91.1 4.638 <.0001
##
## Treatment = MHW2:
## contrast estimate SE df t.ratio p.value
## Black_Head - Palm_Beach 0.12848 0.0362 95.3 3.545 0.0006
##
## Results are averaged over the levels of: Day
## Degrees-of-freedom method: kenward-rogeremmeans(my_lm2, specs = pairwise ~ Day|Treatment, type = "response", adjust = "none")## NOTE: Results may be misleading due to involvement in interactions## $emmeans
## Treatment = Ambient:
## Day emmean SE df lower.CL upper.CL
## 9 0.511 0.0466 26.0 0.415 0.607
## 16 0.713 0.0466 26.0 0.617 0.809
## 23 0.839 0.0466 26.0 0.743 0.934
## 30 0.964 0.0466 26.0 0.868 1.059
##
## Treatment = MHW1:
## Day emmean SE df lower.CL upper.CL
## 9 0.481 0.0466 26.0 0.386 0.577
## 16 0.797 0.0466 26.0 0.701 0.892
## 23 0.971 0.0466 26.0 0.875 1.067
## 30 1.148 0.0466 26.0 1.053 1.244
##
## Treatment = MHW2:
## Day emmean SE df lower.CL upper.CL
## 9 0.507 0.0466 26.0 0.412 0.603
## 16 0.651 0.0519 36.1 0.546 0.757
## 23 0.767 0.0538 40.8 0.658 0.875
## 30 0.878 0.0538 40.8 0.769 0.986
##
## Results are averaged over the levels of: Site
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## Treatment = Ambient:
## contrast estimate SE df t.ratio p.value
## Day9 - Day16 -0.202 0.0393 91.1 -5.144 <.0001
## Day9 - Day23 -0.328 0.0393 91.1 -8.339 <.0001
## Day9 - Day30 -0.453 0.0393 91.1 -11.513 <.0001
## Day16 - Day23 -0.126 0.0393 91.1 -3.195 0.0019
## Day16 - Day30 -0.250 0.0393 91.1 -6.369 <.0001
## Day23 - Day30 -0.125 0.0393 91.1 -3.173 0.0021
##
## Treatment = MHW1:
## contrast estimate SE df t.ratio p.value
## Day9 - Day16 -0.315 0.0393 91.1 -8.015 <.0001
## Day9 - Day23 -0.490 0.0393 91.1 -12.457 <.0001
## Day9 - Day30 -0.667 0.0393 91.1 -16.964 <.0001
## Day16 - Day23 -0.175 0.0393 91.1 -4.441 <.0001
## Day16 - Day30 -0.352 0.0393 91.1 -8.948 <.0001
## Day23 - Day30 -0.177 0.0393 91.1 -4.507 <.0001
##
## Treatment = MHW2:
## contrast estimate SE df t.ratio p.value
## Day9 - Day16 -0.144 0.0455 95.0 -3.169 0.0021
## Day9 - Day23 -0.259 0.0477 94.8 -5.434 <.0001
## Day9 - Day30 -0.370 0.0477 94.8 -7.761 <.0001
## Day16 - Day23 -0.115 0.0508 91.9 -2.266 0.0258
## Day16 - Day30 -0.226 0.0508 91.9 -4.450 <.0001
## Day23 - Day30 -0.111 0.0520 91.1 -2.133 0.0356
##
## Results are averaged over the levels of: Site
## Degrees-of-freedom method: kenward-rogeremmeans(my_lm2, specs = pairwise ~ Treatment|Day, type = "response", adjust = "none")## NOTE: Results may be misleading due to involvement in interactions## $emmeans
## Day = 9:
## Treatment emmean SE df lower.CL upper.CL
## Ambient 0.511 0.0466 26.0 0.415 0.607
## MHW1 0.481 0.0466 26.0 0.386 0.577
## MHW2 0.507 0.0466 26.0 0.412 0.603
##
## Day = 16:
## Treatment emmean SE df lower.CL upper.CL
## Ambient 0.713 0.0466 26.0 0.617 0.809
## MHW1 0.797 0.0466 26.0 0.701 0.892
## MHW2 0.651 0.0519 36.1 0.546 0.757
##
## Day = 23:
## Treatment emmean SE df lower.CL upper.CL
## Ambient 0.839 0.0466 26.0 0.743 0.934
## MHW1 0.971 0.0466 26.0 0.875 1.067
## MHW2 0.767 0.0538 40.8 0.658 0.875
##
## Day = 30:
## Treatment emmean SE df lower.CL upper.CL
## Ambient 0.964 0.0466 26.0 0.868 1.059
## MHW1 1.148 0.0466 26.0 1.053 1.244
## MHW2 0.878 0.0538 40.8 0.769 0.986
##
## Results are averaged over the levels of: Site
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## Day = 9:
## contrast estimate SE df t.ratio p.value
## Ambient - MHW1 0.02947 0.0658 26.0 0.448 0.6582
## Ambient - MHW2 0.00345 0.0658 26.0 0.052 0.9586
## MHW1 - MHW2 -0.02602 0.0658 26.0 -0.395 0.6959
##
## Day = 16:
## contrast estimate SE df t.ratio p.value
## Ambient - MHW1 -0.08344 0.0658 26.0 -1.267 0.2163
## Ambient - MHW2 0.06162 0.0697 31.1 0.884 0.3834
## MHW1 - MHW2 0.14506 0.0697 31.1 2.081 0.0457
##
## Day = 23:
## contrast estimate SE df t.ratio p.value
## Ambient - MHW1 -0.13244 0.0658 26.0 -2.011 0.0548
## Ambient - MHW2 0.07212 0.0712 33.4 1.013 0.3182
## MHW1 - MHW2 0.20456 0.0712 33.4 2.874 0.0070
##
## Day = 30:
## contrast estimate SE df t.ratio p.value
## Ambient - MHW1 -0.18488 0.0658 26.0 -2.808 0.0093
## Ambient - MHW2 0.08593 0.0712 33.4 1.207 0.2358
## MHW1 - MHW2 0.27081 0.0712 33.4 3.805 0.0006
##
## Results are averaged over the levels of: Site
## Degrees-of-freedom method: kenward-roger- Check model assumptions
res <- resid(my_lm2)
#produce residual vs. fitted plot
plot(fitted(my_lm2), res)
#add a horizontal line at 0
abline(0,0)
#create Q-Q plot for residuals
qqnorm(res)
#add a straight diagonal line to the plot
qqline(res) 
2c) Germlings area
1. Visualise differences
Average_Area2 <- Average_Area
Average_Area2$Site <- factor(
Average_Area2$Site,
levels = c("Black_Head", "Palm_Beach"),
labels = c(
"Black Head", "Palm Beach"
))Average_Area2$Day <- as.numeric(Average_Area2$Day)
Area_plot16 <- Average_Area2 %>%
filter(Site == "Palm Beach" | Site == "Black Head") %>%
filter(Day == 9|Day == 16|Day == 23| Day == 30) %>%
ggplot(aes(x= Day, y= MeanArea, colour = Site)) +
geom_errorbar(aes(ymin=MeanArea-se, ymax=MeanArea+se),colour = "grey", width=.2) +
geom_point() +
geom_line()+
Theme1 +
labs( y="Area ("~mm^2*")", x="Time (days)") +
scale_colour_manual(values = c("orange", "darkolivegreen3")) +
theme(strip.text.x = element_text(size = 12))+
scale_y_continuous(limits = c(0, 0.4), breaks = seq(0, 0.4, by = 0.1))+
scale_x_continuous(limits = c(8, 31), breaks = c(9, 16, 23,30))+
facet_wrap(~ factor(Treatment, c("Ambient","MHW1","MHW2")), ncol = 3)
Area_plot16
ggsave("Area_germlings.jpeg", width = 200, height = 56, units = "mm")2. Statistical analyses
- Choose only the days to analyse and transform them to factor
Area_well1 <- filter (Area_tank, Site == "Palm_Beach" | Site == "Black_Head")
Area_well1$Day <- as.factor(Area_well1$Day)
Area_wellALL <- filter (Area_well1, Day == "9" |Day == "16"|Day == "23" |Day == "30")- Run linear model with site, MHW treatment and day as fixes/crossed factors and with water bath (Tank) as random intercept.
my_lm2 <- lmer(Well_meanArea ~ Site *Treatment*Day + (1|TankCode), data=Area_wellALL)
my_lm_anova2 <- Anova(my_lm2)
my_lm_anova2## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Well_meanArea
## Chisq Df Pr(>Chisq)
## Site 29.6580 1 5.154e-08 ***
## Treatment 12.0648 2 0.002400 **
## Day 882.6266 3 < 2.2e-16 ***
## Site:Treatment 3.0722 2 0.215222
## Site:Day 13.7966 3 0.003195 **
## Treatment:Day 50.9940 6 2.970e-09 ***
## Site:Treatment:Day 8.5246 6 0.202129
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- Post-hoc test
emmeans(my_lm2, specs = pairwise ~ Treatment|Day, type = "response", adjust = "none")## NOTE: Results may be misleading due to involvement in interactions## $emmeans
## Day = 9:
## Treatment emmean SE df lower.CL upper.CL
## Ambient 0.0728 0.0115 29.1 0.0494 0.0963
## MHW1 0.0652 0.0115 29.1 0.0417 0.0886
## MHW2 0.0731 0.0115 29.1 0.0497 0.0966
##
## Day = 16:
## Treatment emmean SE df lower.CL upper.CL
## Ambient 0.1394 0.0115 29.1 0.1160 0.1628
## MHW1 0.1646 0.0115 29.1 0.1411 0.1880
## MHW2 0.1149 0.0130 41.2 0.0887 0.1411
##
## Day = 23:
## Treatment emmean SE df lower.CL upper.CL
## Ambient 0.1880 0.0115 29.1 0.1645 0.2114
## MHW1 0.2290 0.0115 29.1 0.2056 0.2525
## MHW2 0.1489 0.0135 46.9 0.1217 0.1762
##
## Day = 30:
## Treatment emmean SE df lower.CL upper.CL
## Ambient 0.2387 0.0115 29.1 0.2153 0.2621
## MHW1 0.2998 0.0115 29.1 0.2764 0.3233
## MHW2 0.2042 0.0135 46.9 0.1770 0.2314
##
## Results are averaged over the levels of: Site
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## Day = 9:
## contrast estimate SE df t.ratio p.value
## Ambient - MHW1 0.007662 0.0162 29.1 0.473 0.6401
## Ambient - MHW2 -0.000288 0.0162 29.1 -0.018 0.9859
## MHW1 - MHW2 -0.007950 0.0162 29.1 -0.490 0.6276
##
## Day = 16:
## contrast estimate SE df t.ratio p.value
## Ambient - MHW1 -0.025159 0.0162 29.1 -1.552 0.1316
## Ambient - MHW2 0.024517 0.0173 35.2 1.417 0.1654
## MHW1 - MHW2 0.049676 0.0173 35.2 2.870 0.0069
##
## Day = 23:
## contrast estimate SE df t.ratio p.value
## Ambient - MHW1 -0.041040 0.0162 29.1 -2.531 0.0170
## Ambient - MHW2 0.039017 0.0177 38.1 2.201 0.0339
## MHW1 - MHW2 0.080057 0.0177 38.1 4.516 0.0001
##
## Day = 30:
## contrast estimate SE df t.ratio p.value
## Ambient - MHW1 -0.061107 0.0162 29.1 -3.769 0.0007
## Ambient - MHW2 0.034497 0.0177 38.1 1.946 0.0591
## MHW1 - MHW2 0.095604 0.0177 38.1 5.393 <.0001
##
## Results are averaged over the levels of: Site
## Degrees-of-freedom method: kenward-roger emmeans(my_lm2, specs = pairwise ~ Day|Treatment, type = "response", adjust = "none")## NOTE: Results may be misleading due to involvement in interactions## $emmeans
## Treatment = Ambient:
## Day emmean SE df lower.CL upper.CL
## 9 0.0728 0.0115 29.1 0.0494 0.0963
## 16 0.1394 0.0115 29.1 0.1160 0.1628
## 23 0.1880 0.0115 29.1 0.1645 0.2114
## 30 0.2387 0.0115 29.1 0.2153 0.2621
##
## Treatment = MHW1:
## Day emmean SE df lower.CL upper.CL
## 9 0.0652 0.0115 29.1 0.0417 0.0886
## 16 0.1646 0.0115 29.1 0.1411 0.1880
## 23 0.2290 0.0115 29.1 0.2056 0.2525
## 30 0.2998 0.0115 29.1 0.2764 0.3233
##
## Treatment = MHW2:
## Day emmean SE df lower.CL upper.CL
## 9 0.0731 0.0115 29.1 0.0497 0.0966
## 16 0.1149 0.0130 41.2 0.0887 0.1411
## 23 0.1489 0.0135 46.9 0.1217 0.1762
## 30 0.2042 0.0135 46.9 0.1770 0.2314
##
## Results are averaged over the levels of: Site
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## Treatment = Ambient:
## contrast estimate SE df t.ratio p.value
## Day9 - Day16 -0.0666 0.0105 91.1 -6.332 <.0001
## Day9 - Day23 -0.1151 0.0105 91.1 -10.953 <.0001
## Day9 - Day30 -0.1659 0.0105 91.1 -15.780 <.0001
## Day16 - Day23 -0.0486 0.0105 91.1 -4.621 <.0001
## Day16 - Day30 -0.0993 0.0105 91.1 -9.448 <.0001
## Day23 - Day30 -0.0507 0.0105 91.1 -4.827 <.0001
##
## Treatment = MHW1:
## contrast estimate SE df t.ratio p.value
## Day9 - Day16 -0.0994 0.0105 91.1 -9.455 <.0001
## Day9 - Day23 -0.1638 0.0105 91.1 -15.587 <.0001
## Day9 - Day30 -0.2346 0.0105 91.1 -22.323 <.0001
## Day16 - Day23 -0.0644 0.0105 91.1 -6.132 <.0001
## Day16 - Day30 -0.1352 0.0105 91.1 -12.868 <.0001
## Day23 - Day30 -0.0708 0.0105 91.1 -6.736 <.0001
##
## Treatment = MHW2:
## contrast estimate SE df t.ratio p.value
## Day9 - Day16 -0.0417 0.0121 95.5 -3.443 0.0009
## Day9 - Day23 -0.0758 0.0127 95.3 -5.959 <.0001
## Day9 - Day30 -0.1311 0.0127 95.3 -10.303 <.0001
## Day16 - Day23 -0.0341 0.0136 92.1 -2.509 0.0139
## Day16 - Day30 -0.0893 0.0136 92.1 -6.579 <.0001
## Day23 - Day30 -0.0553 0.0139 91.1 -3.974 0.0001
##
## Results are averaged over the levels of: Site
## Degrees-of-freedom method: kenward-rogeremmeans(my_lm2, specs = pairwise ~ Site|Day, type = "response", adjust = "none")## NOTE: Results may be misleading due to involvement in interactions## $emmeans
## Day = 9:
## Site emmean SE df lower.CL upper.CL
## Black_Head 0.0714 0.00789 51.7 0.0555 0.0872
## Palm_Beach 0.0694 0.00789 51.7 0.0536 0.0852
##
## Day = 16:
## Site emmean SE df lower.CL upper.CL
## Black_Head 0.1486 0.00835 59.0 0.1319 0.1653
## Palm_Beach 0.1307 0.00835 58.9 0.1139 0.1474
##
## Day = 23:
## Site emmean SE df lower.CL upper.CL
## Black_Head 0.2064 0.00835 59.0 0.1897 0.2231
## Palm_Beach 0.1709 0.00877 65.7 0.1534 0.1884
##
## Day = 30:
## Site emmean SE df lower.CL upper.CL
## Black_Head 0.2733 0.00835 59.0 0.2566 0.2900
## Palm_Beach 0.2218 0.00877 65.7 0.2043 0.2394
##
## Results are averaged over the levels of: Treatment
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## Day = 9:
## contrast estimate SE df t.ratio p.value
## Black_Head - Palm_Beach 0.00196 0.00858 91.1 0.228 0.8198
##
## Day = 16:
## contrast estimate SE df t.ratio p.value
## Black_Head - Palm_Beach 0.01792 0.00935 91.8 1.917 0.0583
##
## Day = 23:
## contrast estimate SE df t.ratio p.value
## Black_Head - Palm_Beach 0.03552 0.00976 92.7 3.640 0.0004
##
## Day = 30:
## contrast estimate SE df t.ratio p.value
## Black_Head - Palm_Beach 0.05144 0.00976 92.7 5.272 <.0001
##
## Results are averaged over the levels of: Treatment
## Degrees-of-freedom method: kenward-rogeremmeans(my_lm2, specs = pairwise ~ Day|Site, type = "response", adjust = "none")## NOTE: Results may be misleading due to involvement in interactions## $emmeans
## Site = Black_Head:
## Day emmean SE df lower.CL upper.CL
## 9 0.0714 0.00789 51.7 0.0555 0.0872
## 16 0.1486 0.00835 59.0 0.1319 0.1653
## 23 0.2064 0.00835 59.0 0.1897 0.2231
## 30 0.2733 0.00835 59.0 0.2566 0.2900
##
## Site = Palm_Beach:
## Day emmean SE df lower.CL upper.CL
## 9 0.0694 0.00789 51.7 0.0536 0.0852
## 16 0.1307 0.00835 58.9 0.1139 0.1474
## 23 0.1709 0.00877 65.7 0.1534 0.1884
## 30 0.2218 0.00877 65.7 0.2043 0.2394
##
## Results are averaged over the levels of: Treatment
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## Site = Black_Head:
## contrast estimate SE df t.ratio p.value
## Day9 - Day16 -0.0772 0.00900 92.2 -8.575 <.0001
## Day9 - Day23 -0.1350 0.00901 92.2 -14.990 <.0001
## Day9 - Day30 -0.2019 0.00901 92.2 -22.415 <.0001
## Day16 - Day23 -0.0578 0.00934 91.7 -6.192 <.0001
## Day16 - Day30 -0.1247 0.00934 91.7 -13.353 <.0001
## Day23 - Day30 -0.0669 0.00927 91.1 -7.216 <.0001
##
## Site = Palm_Beach:
## contrast estimate SE df t.ratio p.value
## Day9 - Day16 -0.0612 0.00901 92.3 -6.796 <.0001
## Day9 - Day23 -0.1015 0.00940 92.7 -10.796 <.0001
## Day9 - Day30 -0.1524 0.00940 92.7 -16.218 <.0001
## Day16 - Day23 -0.0402 0.00964 91.5 -4.172 0.0001
## Day16 - Day30 -0.0912 0.00964 91.5 -9.457 <.0001
## Day23 - Day30 -0.0510 0.00991 91.1 -5.144 <.0001
##
## Results are averaged over the levels of: Treatment
## Degrees-of-freedom method: kenward-roger- Check model assumptions
res <- resid(my_lm2)
#produce residual vs. fitted plot
plot(fitted(my_lm2), res)
#add a horizontal line at 0
abline(0,0)
#create Q-Q plot for residuals
qqnorm(res)
#add a straight diagonal line to the plot
qqline(res) 
2d) Shape
1. Visualise differences
Average_Shape2 <- Average_Shape
Average_Shape2$Site <- factor(
Average_Shape$Site,
levels = c("Palm_Beach", "Black_Head"),
labels = c(
"Palm Beach",
"Black Head"
))Average_Shape2 <- Average_Shape2 %>%
filter(Site == "Palm Beach" | Site == "Black Head") %>%
filter(Day == 9|Day == 16|Day == 23| Day == 30)
Average_Shape2$Day <- as.numeric(Average_Shape2$Day)
Shapes <- Average_Shape2 %>%
ggplot(aes(x= Day, y= MeanShape, colour = Site)) +
geom_errorbar(aes(ymin=MeanShape-se, ymax=MeanShape+se),colour = "grey", width=.2) +
geom_point() +
geom_line()+
Theme1 +
labs( y="Proportion of occurrence", x="Time (days)") +
scale_colour_manual(values = c("orange", "darkolivegreen3")) +
theme(strip.text.x = element_text(size = 12))+
scale_y_continuous(limits = c(0, 1), breaks = seq(0, 1, by = 0.2))+
scale_x_continuous(limits = c(8, 31), breaks = c(9, 16, 23,30))+
facet_wrap(~ Shape2 +factor(Treatment, c("Ambient","MHW1","MHW2")), ncol = 3)
Shapes
ggsave("shape_germlings.jpeg", width = 180, height = 140, units = "mm")2. Statistical analyses
Separate by shape and select the days to analyse
Shape_long <- filter (Shape_long, Site == "Palm_Beach" | Site == "Black_Head")
Shape_long$Day <- as.factor(Shape_long$Day)
Shape_el <- filter (Shape_long, Shape2 == "Elongated")
Shape_well1 <- filter (Shape_el, Day == "9" |Day == "16"|Day == "23" |Day == "30")
Shape_div <- filter (Shape_long, Shape2 == "Divided")
ShapeD_well1 <- filter (Shape_div, Day == "9" |Day == "16"|Day == "23" |Day == "30")
Shape_r <- filter (Shape_long, Shape2 == "Rounded")
ShapeR_well1 <- filter (Shape_r, Day == "9" |Day == "16"|Day == "23" |Day == "30")2.1. Branched (divided)
- Run LMM for percentage of branched germlings
my_lm2 <- lmer(percent ~ Site *Treatment*Day + (1|TankCode), data=ShapeD_well1)
summary(my_lm2)## Linear mixed model fit by REML ['lmerMod']
## Formula: percent ~ Site * Treatment * Day + (1 | TankCode)
## Data: ShapeD_well1
##
## REML criterion at convergence: 28.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.1997 -0.4598 -0.1989 0.4701 2.7282
##
## Random effects:
## Groups Name Variance Std.Dev.
## TankCode (Intercept) 0.0007258 0.02694
## Residual 0.0518397 0.22768
## Number of obs: 130, groups: TankCode, 18
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 3.645e-15 9.360e-02 0.000
## SitePalm_Beach 1.726e-01 1.315e-01 1.313
## TreatmentMHW1 5.556e-02 1.324e-01 0.420
## TreatmentMHW2 5.556e-02 1.324e-01 0.420
## Day16 1.087e-01 1.315e-01 0.827
## Day23 3.250e-01 1.315e-01 2.472
## Day30 3.111e-01 1.315e-01 2.367
## SitePalm_Beach:TreatmentMHW1 -1.671e-01 1.859e-01 -0.899
## SitePalm_Beach:TreatmentMHW2 -1.726e-01 1.859e-01 -0.929
## SitePalm_Beach:Day16 -4.524e-02 1.859e-01 -0.243
## SitePalm_Beach:Day23 -1.683e-01 1.859e-01 -0.905
## SitePalm_Beach:Day30 4.881e-02 1.859e-01 0.263
## TreatmentMHW1:Day16 -3.929e-02 1.859e-01 -0.211
## TreatmentMHW2:Day16 -7.931e-02 1.973e-01 -0.402
## TreatmentMHW1:Day23 -2.778e-03 1.859e-01 -0.015
## TreatmentMHW2:Day23 -2.993e-01 1.973e-01 -1.517
## TreatmentMHW1:Day30 7.778e-02 1.859e-01 0.418
## TreatmentMHW2:Day30 1.313e-01 1.973e-01 0.665
## SitePalm_Beach:TreatmentMHW1:Day16 1.647e-01 2.629e-01 0.626
## SitePalm_Beach:TreatmentMHW2:Day16 1.300e-01 2.790e-01 0.466
## SitePalm_Beach:TreatmentMHW1:Day23 2.036e-01 2.629e-01 0.774
## SitePalm_Beach:TreatmentMHW2:Day23 3.013e-01 2.868e-01 1.051
## SitePalm_Beach:TreatmentMHW1:Day30 -6.548e-02 2.629e-01 -0.249
## SitePalm_Beach:TreatmentMHW2:Day30 -4.356e-01 2.868e-01 -1.519##
## Correlation matrix not shown by default, as p = 24 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need itAnova(my_lm2)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: percent
## Chisq Df Pr(>Chisq)
## Site 1.9257 1 0.1652
## Treatment 3.4810 2 0.1754
## Day 41.1430 3 6.098e-09 ***
## Site:Treatment 2.7048 2 0.2586
## Site:Day 1.2014 3 0.7527
## Treatment:Day 3.5466 6 0.7378
## Site:Treatment:Day 6.5824 6 0.3612
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- Post-hoc test
emm5 = emmeans(my_lm2, specs = pairwise ~ Day, type = "response", adjust = "none")## NOTE: Results may be misleading due to involvement in interactionsemm5## $emmeans
## Day emmean SE df lower.CL upper.CL
## 9 0.0667 0.0385 88.0 -0.00973 0.143
## 16 0.1624 0.0416 91.3 0.07980 0.245
## 23 0.2911 0.0430 93.7 0.20559 0.377
## 30 0.3884 0.0430 93.7 0.30295 0.474
##
## Results are averaged over the levels of: Site, Treatment
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## Day9 - Day16 -0.0957 0.0559 94.1 -1.711 0.0905
## Day9 - Day23 -0.2243 0.0570 94.4 -3.933 0.0002
## Day9 - Day30 -0.3217 0.0570 94.4 -5.640 <.0001
## Day16 - Day23 -0.1286 0.0591 93.7 -2.176 0.0321
## Day16 - Day30 -0.2260 0.0591 93.7 -3.823 0.0002
## Day23 - Day30 -0.0974 0.0600 91.8 -1.623 0.1081
##
## Results are averaged over the levels of: Site, Treatment
## Degrees-of-freedom method: kenward-roger- Check model assumptions
res <- resid(my_lm2)
#produce residual vs. fitted plot
plot(fitted(my_lm2), res)
#add a horizontal line at 0
abline(0,0)
#create Q-Q plot for residuals
qqnorm(res)
#add a straight diagonal line to the plot
qqline(res) 
2.2. Elongated
- Run LMM for percentage of elongated germlings
my_lm2 <- lmer(percent ~ Site *Treatment*Day + (1|TankCode), data=Shape_well1)
summary(my_lm2)## Linear mixed model fit by REML ['lmerMod']
## Formula: percent ~ Site * Treatment * Day + (1 | TankCode)
## Data: Shape_well1
##
## REML criterion at convergence: 65.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.09880 -0.59045 -0.01315 0.62542 1.97960
##
## Random effects:
## Groups Name Variance Std.Dev.
## TankCode (Intercept) 0.008365 0.09146
## Residual 0.068313 0.26137
## Number of obs: 130, groups: TankCode, 18
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.77778 0.11305 6.880
## SitePalm_Beach -0.07817 0.15090 -0.518
## TreatmentMHW1 -0.18333 0.15987 -1.147
## TreatmentMHW2 -0.10833 0.15987 -0.678
## Day16 -0.10476 0.15090 -0.694
## Day23 -0.22778 0.15090 -1.509
## Day30 -0.27500 0.15090 -1.822
## SitePalm_Beach:TreatmentMHW1 0.05040 0.21341 0.236
## SitePalm_Beach:TreatmentMHW2 0.14206 0.21341 0.666
## SitePalm_Beach:Day16 0.07183 0.21341 0.337
## SitePalm_Beach:Day23 0.19881 0.21341 0.932
## SitePalm_Beach:Day30 0.04286 0.21341 0.201
## TreatmentMHW1:Day16 0.21865 0.21341 1.025
## TreatmentMHW2:Day16 -0.05253 0.22736 -0.231
## TreatmentMHW1:Day23 0.11667 0.21341 0.547
## TreatmentMHW2:Day23 0.39036 0.22740 1.717
## TreatmentMHW1:Day30 0.18056 0.21341 0.846
## TreatmentMHW2:Day30 0.10425 0.22740 0.458
## SitePalm_Beach:TreatmentMHW1:Day16 -0.11349 0.30180 -0.376
## SitePalm_Beach:TreatmentMHW2:Day16 0.03812 0.32109 0.119
## SitePalm_Beach:TreatmentMHW1:Day23 -0.10079 0.30180 -0.334
## SitePalm_Beach:TreatmentMHW2:Day23 -0.37345 0.33087 -1.129
## SitePalm_Beach:TreatmentMHW1:Day30 0.05159 0.30180 0.171
## SitePalm_Beach:TreatmentMHW2:Day30 -0.06670 0.33087 -0.202##
## Correlation matrix not shown by default, as p = 24 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need itAnova(my_lm2)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: percent
## Chisq Df Pr(>Chisq)
## Site 0.1305 1 0.7179
## Treatment 0.8582 2 0.6511
## Day 7.0093 3 0.0716 .
## Site:Treatment 0.2492 2 0.8829
## Site:Day 0.3061 3 0.9589
## Treatment:Day 6.1533 6 0.4062
## Site:Treatment:Day 2.3889 6 0.8807
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- Post-hoc test
emm5 = emmeans(my_lm2, specs = pairwise ~ Day, type = "response", adjust = "none")## NOTE: Results may be misleading due to involvement in interactionsemm5## $emmeans
## Day emmean SE df lower.CL upper.CL
## 9 0.674 0.0486 63.4 0.576 0.771
## 16 0.648 0.0522 69.7 0.543 0.752
## 23 0.635 0.0537 73.5 0.528 0.742
## 30 0.512 0.0537 73.5 0.405 0.619
##
## Results are averaged over the levels of: Site, Treatment
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## Day9 - Day16 0.0260 0.0645 93.8 0.404 0.6872
## Day9 - Day23 0.0384 0.0657 94.0 0.584 0.5605
## Day9 - Day30 0.1612 0.0657 94.0 2.452 0.0161
## Day16 - Day23 0.0124 0.0680 92.8 0.182 0.8561
## Day16 - Day30 0.1351 0.0680 92.8 1.987 0.0498
## Day23 - Day30 0.1228 0.0689 91.5 1.782 0.0780
##
## Results are averaged over the levels of: Site, Treatment
## Degrees-of-freedom method: kenward-roger- Check model assumptions
res <- resid(my_lm2)
#produce residual vs. fitted plot
plot(fitted(my_lm2), res)
#add a horizontal line at 0
abline(0,0)
#create Q-Q plot for residuals
qqnorm(res)
#add a straight diagonal line to the plot
qqline(res) 
2.3. Rounded
- Run LMM for percentage of rounded germlings
my_lm2 <- lmer(percent ~ Site *Treatment*Day + (1|TankCode), data=ShapeR_well1)
summary(my_lm2)## Linear mixed model fit by REML ['lmerMod']
## Formula: percent ~ Site * Treatment * Day + (1 | TankCode)
## Data: ShapeR_well1
##
## REML criterion at convergence: -4.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8955 -0.4793 -0.0982 0.4476 3.2300
##
## Random effects:
## Groups Name Variance Std.Dev.
## TankCode (Intercept) 0.01209 0.1099
## Residual 0.03205 0.1790
## Number of obs: 130, groups: TankCode, 18
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.222222 0.085774 2.591
## SitePalm_Beach -0.094444 0.103367 -0.914
## TreatmentMHW1 0.127778 0.121303 1.053
## TreatmentMHW2 0.052778 0.121303 0.435
## Day16 -0.003968 0.103367 -0.038
## Day23 -0.097222 0.103367 -0.941
## Day30 -0.036111 0.103367 -0.349
## SitePalm_Beach:TreatmentMHW1 0.116667 0.146183 0.798
## SitePalm_Beach:TreatmentMHW2 0.030556 0.146183 0.209
## SitePalm_Beach:Day16 -0.026587 0.146183 -0.182
## SitePalm_Beach:Day23 -0.030556 0.146183 -0.209
## SitePalm_Beach:Day30 -0.091667 0.146183 -0.627
## TreatmentMHW1:Day16 -0.179365 0.146183 -1.227
## TreatmentMHW2:Day16 0.124923 0.156356 0.799
## TreatmentMHW1:Day23 -0.113889 0.146183 -0.779
## TreatmentMHW2:Day23 -0.073130 0.156422 -0.468
## TreatmentMHW1:Day30 -0.258333 0.146183 -1.767
## TreatmentMHW2:Day30 -0.217575 0.156422 -1.391
## SitePalm_Beach:TreatmentMHW1:Day16 -0.051190 0.206734 -0.248
## SitePalm_Beach:TreatmentMHW2:Day16 -0.173338 0.220419 -0.786
## SitePalm_Beach:TreatmentMHW1:Day23 -0.102778 0.206734 -0.497
## SitePalm_Beach:TreatmentMHW2:Day23 0.061443 0.227643 0.270
## SitePalm_Beach:TreatmentMHW1:Day30 0.013889 0.206734 0.067
## SitePalm_Beach:TreatmentMHW2:Day30 0.491602 0.227643 2.160##
## Correlation matrix not shown by default, as p = 24 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need itAnova(my_lm2)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: percent
## Chisq Df Pr(>Chisq)
## Site 5.6784 1 0.01717 *
## Treatment 0.5534 2 0.75826
## Day 22.1363 3 6.111e-05 ***
## Site:Treatment 1.8918 2 0.38834
## Site:Day 1.8542 3 0.60321
## Treatment:Day 9.5694 6 0.14400
## Site:Treatment:Day 9.8251 6 0.13221
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- Post-hoc test
emm5 = emmeans(my_lm2, specs = pairwise ~ Day, type = "response", adjust = "none")## NOTE: Results may be misleading due to involvement in interactionsemm5## $emmeans
## Day emmean SE df lower.CL upper.CL
## 9 0.260 0.0395 38.9 0.17977 0.340
## 16 0.187 0.0418 44.6 0.10265 0.271
## 23 0.078 0.0427 47.5 -0.00798 0.164
## 30 0.103 0.0427 47.5 0.01741 0.189
##
## Results are averaged over the levels of: Site, Treatment
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## Day9 - Day16 0.0728 0.0444 93.3 1.642 0.1040
## Day9 - Day23 0.1817 0.0452 93.4 4.017 0.0001
## Day9 - Day30 0.1563 0.0452 93.4 3.456 0.0008
## Day16 - Day23 0.1089 0.0466 91.9 2.335 0.0217
## Day16 - Day30 0.0835 0.0466 91.9 1.790 0.0767
## Day23 - Day30 -0.0254 0.0472 91.2 -0.538 0.5917
##
## Results are averaged over the levels of: Site, Treatment
## Degrees-of-freedom method: kenward-rogeremm5 = emmeans(my_lm2, specs = pairwise ~ Site, type = "response", adjust = "none")## NOTE: Results may be misleading due to involvement in interactionsemm5## $emmeans
## Site emmean SE df lower.CL upper.CL
## Black_Head 0.188 0.0346 23.4 0.1168 0.260
## Palm_Beach 0.126 0.0354 25.0 0.0528 0.199
##
## Results are averaged over the levels of: Treatment, Day
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## Black_Head - Palm_Beach 0.0626 0.033 94.7 1.897 0.0608
##
## Results are averaged over the levels of: Treatment, Day
## Degrees-of-freedom method: kenward-roger- Check model assumptions
res <- resid(my_lm2)
#produce residual vs. fitted plot
plot(fitted(my_lm2), res)
#add a horizontal line at 0
abline(0,0)
#create Q-Q plot for residuals
qqnorm(res)
#add a straight diagonal line to the plot
qqline(res) 
Section 3: Juvenile experiments
3a) Survival
1. Visualise differences
surv_plot <- Juv_means %>%
filter(Day_equivalent == 8|Day_equivalent == 12|Day_equivalent == 18| Day_equivalent == 22) %>%
ggplot(aes(x= Day_equivalent, y= Mean_survival, colour = Site)) +
geom_errorbar(aes(ymin=Mean_survival-se_surv, ymax=Mean_survival+se_surv),colour = "grey", width=.2) +
geom_point() +
geom_line()+
Theme1 +
labs( y="Survival (0 or 1)", x="Time (days)") +
scale_colour_manual(values = c("darkred", "orange", "darkolivegreen3", "royalblue3")) +
theme(strip.text.x = element_text(size = 12))+
scale_y_continuous(limits = c(0, 1.05), breaks = seq(0, 1, by = 1))+
scale_x_continuous(limits = c(7, 23), breaks = c(8, 12, 18,22))+
facet_wrap(~ factor(Treatment, c("Ambient","MHW1","MHW2")), ncol = 3)
surv_plot
ggsave("Survival01_juv.jpeg", width = 220, height = 60, units = "mm")2. Statistical analyses
- Separate data by day
Juveniles_meta$Tank_Code <- as.factor(Juveniles_meta$Tank_Code)
Juveniles_meta$Day_equivalent <- as.factor(Juveniles_meta$Day_equivalent)
##Filter Juveniles_meta for survival and holdfast analyses
Juveniles_stat8_1<- Juveniles_meta%>%
filter( Day_equivalent == "8")
Juveniles_stat12_1<- Juveniles_meta%>%
filter( Day_equivalent == "12")
Juveniles_stat18_1<- Juveniles_meta%>%
filter( Day_equivalent == "18")
Juveniles_stat22_1<- Juveniles_meta%>%
filter( Day_equivalent == "22")
Juveniles_stat1<- Juveniles_meta %>%
filter(Day_equivalent == "8"|Day_equivalent == "12"| Day_equivalent == "18"| Day_equivalent == "22")2.1. Day 12
- Run LMM for survival at day 12
my_glmer<- glmer(Survival_code ~ Site*Treatment + (1|Location), data=Juveniles_stat12_1, family = binomial)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## unable to evaluate scaled gradient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Hessian is numerically singular: parameters are not uniquely determinedsummary(my_glmer)## Warning in vcov.merMod(object, use.hessian = use.hessian): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX## Warning in vcov.merMod(object, correlation = correlation, sigm = sig): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Survival_code ~ Site * Treatment + (1 | Location)
## Data: Juveniles_stat12_1
##
## AIC BIC logLik deviance df.resid
## 51.6 88.0 -12.8 25.6 109
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.92019 0.00001 0.00002 0.24255 0.45208
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 0.4635 0.6808
## Number of obs: 122, groups: Location, 6
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.273e+01 2.977e+04 0.001 0.999
## SiteBonny Hills -1.980e-01 3.973e+04 0.000 1.000
## SiteCronulla -2.014e+01 2.977e+04 -0.001 0.999
## SitePalm Beach 5.389e-01 4.436e+04 0.000 1.000
## TreatmentMHW1 5.481e+00 3.838e+05 0.000 1.000
## TreatmentMHW2 -2.041e+01 2.977e+04 -0.001 0.999
## SiteBonny Hills:TreatmentMHW1 -1.822e+00 4.195e+05 0.000 1.000
## SiteCronulla:TreatmentMHW1 1.321e+01 3.840e+05 0.000 1.000
## SitePalm Beach:TreatmentMHW1 -6.823e+00 3.856e+05 0.000 1.000
## SiteBonny Hills:TreatmentMHW2 -5.943e-02 3.973e+04 0.000 1.000
## SiteCronulla:TreatmentMHW2 3.855e+01 3.104e+04 0.001 0.999
## SitePalm Beach:TreatmentMHW2 -6.994e-01 4.436e+04 0.000 1.000
##
## Correlation of Fixed Effects:
## (Intr) StBnnH StCrnl StPlmB TrMHW1 TrMHW2 SBH:TMHW1 SC:TMHW1
## SitBnnyHlls -0.749
## SiteCronull -1.000 0.749
## SitePalmBch -0.671 0.503 0.671
## TretmntMHW1 -0.078 0.058 0.078 0.052
## TretmntMHW2 -1.000 0.749 1.000 0.671 0.078
## StBHl:TMHW1 0.071 -0.095 -0.071 -0.048 -0.915 -0.071
## StCrn:TMHW1 0.078 -0.058 -0.078 -0.052 -1.000 -0.078 0.914
## StPBc:TMHW1 0.077 -0.058 -0.077 -0.115 -0.995 -0.077 0.911 0.995
## StBHl:TMHW2 0.749 -1.000 -0.749 -0.503 -0.058 -0.749 0.095 0.058
## StCrn:TMHW2 0.959 -0.718 -0.959 -0.643 -0.074 -0.959 0.068 0.074
## StPBc:TMHW2 0.671 -0.503 -0.671 -1.000 -0.052 -0.671 0.048 0.052
## SPB:TMHW1 SBH:TMHW2 SC:TMHW2
## SitBnnyHlls
## SiteCronull
## SitePalmBch
## TretmntMHW1
## TretmntMHW2
## StBHl:TMHW1
## StCrn:TMHW1
## StPBc:TMHW1
## StBHl:TMHW2 0.058
## StCrn:TMHW2 0.074 0.718
## StPBc:TMHW2 0.115 0.503 0.643
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## unable to evaluate scaled gradient
## Hessian is numerically singular: parameters are not uniquely determinedmy_glmer_anova <- car::Anova(my_glmer)## Warning in vcov.merMod(mod, complete = FALSE): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RXmy_glmer_anova## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Survival_code
## Chisq Df Pr(>Chisq)
## Site 0.0274 3 0.9988
## Treatment 0.0000 2 1.0000
## Site:Treatment 0.0000 6 1.0000overdisp.glmer(my_glmer)## Residual deviance: 23.031 on 109 degrees of freedom (ratio: 0.211)- Check model assumptions
res <- resid(my_glmer)
#produce residual vs. fitted plot
plot(fitted(my_glmer), res)
#add a horizontal line at 0
abline(0,0)
#create Q-Q plot for residuals
qqnorm(res)
#add a straight diagonal line to the plot
qqline(res) 
2.2. Day 18
- Run LMM for survival at day 18
my_glmer1<- glmer(Survival_code ~ Site*Treatment + (1|Location), data=Juveniles_stat18_1, family = binomial)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Hessian is numerically singular: parameters are not uniquely determinedsummary(my_glmer1)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Survival_code ~ Site * Treatment + (1 | Location)
## Data: Juveniles_stat18_1
##
## AIC BIC logLik deviance df.resid
## 84.8 121.3 -29.4 58.8 109
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.07651 -0.00003 0.00003 0.40997 2.33026
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 0.3544 0.5953
## Number of obs: 122, groups: Location, 6
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.2346 1.1266 1.984 0.0473 *
## SiteBonny Hills 18.3283 340.4831 0.054 0.9571
## SiteCronulla -0.5097 1.3319 -0.383 0.7019
## SitePalm Beach -0.6213 1.3394 -0.464 0.6427
## TreatmentMHW1 -0.6117 1.3386 -0.457 0.6477
## TreatmentMHW2 -22.9213 846.9958 -0.027 0.9784
## SiteBonny Hills:TreatmentMHW1 -18.8558 340.4837 -0.055 0.9558
## SiteCronulla:TreatmentMHW1 18.6373 628.5206 0.030 0.9763
## SitePalm Beach:TreatmentMHW1 20.6494 519.1645 0.040 0.9683
## SiteBonny Hills:TreatmentMHW2 0.1482 808.5563 0.000 0.9999
## SiteCronulla:TreatmentMHW2 -0.4994 857.7094 -0.001 0.9995
## SitePalm Beach:TreatmentMHW2 19.0934 846.9961 0.023 0.9820
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) StBnnH StCrnl StPlmB TrMHW1 TrMHW2 SBH:TMHW1 SC:TMHW1
## SitBnnyHlls -0.001
## SiteCronull -0.772 0.001
## SitePalmBch -0.770 0.001 0.646
## TretmntMHW1 -0.768 0.000 0.646 0.643
## TretmntMHW2 0.001 -0.311 -0.001 -0.002 -0.001
## StBHl:TMHW1 0.001 -1.000 -0.001 -0.001 -0.002 0.311
## StCrn:TMHW1 -0.001 -0.138 0.001 0.001 0.001 -0.208 0.138
## StPBc:TMHW1 0.000 0.157 0.000 0.000 0.000 -0.074 -0.157 -0.044
## StBHl:TMHW2 -0.002 -0.095 0.002 0.003 0.002 -0.916 0.095 0.276
## StCrn:TMHW2 0.002 -0.155 -0.001 -0.002 -0.002 0.711 0.155 -0.274
## StPBc:TMHW2 -0.002 0.311 0.001 0.001 0.001 -1.000 -0.311 0.208
## SPB:TMHW1 SBH:TMHW2 SC:TMHW2
## SitBnnyHlls
## SiteCronull
## SitePalmBch
## TretmntMHW1
## TretmntMHW2
## StBHl:TMHW1
## StCrn:TMHW1
## StPBc:TMHW1
## StBHl:TMHW2 0.011
## StCrn:TMHW2 -0.064 -0.680
## StPBc:TMHW2 0.074 0.916 -0.711
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## Hessian is numerically singular: parameters are not uniquely determinedoverdisp.glmer(my_glmer1)## Residual deviance: 54.354 on 109 degrees of freedom (ratio: 0.499)my_glmer_anova <- car::Anova(my_glmer1)
my_glmer_anova## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Survival_code
## Chisq Df Pr(>Chisq)
## Site 0.4211 3 0.935858
## Treatment 11.7712 2 0.002779 **
## Site:Treatment 0.0076 6 1.000000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- Post-hoc test
emm3 <- emmeans(my_glmer1, specs = pairwise ~ Treatment, type = "response", adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionssummary(emm3)## $emmeans
## Treatment prob SE df asymp.LCL asymp.UCL
## Ambient 9.99e-01 0.12336 Inf 2.22e-16 1
## MHW1 1.00e+00 0.00323 Inf 2.22e-16 1
## MHW2 8.30e-06 0.00493 Inf 2.22e-16 1
##
## Results are averaged over the levels of: Site
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scale
##
## $contrasts
## contrast odds.ratio SE df null z.ratio p.value
## Ambient / MHW1 0.00e+00 2.00e+00 Inf 1 -0.021 0.9998
## Ambient / MHW2 8.31e+07 5.19e+10 Inf 1 0.029 0.9995
## MHW1 / MHW2 7.45e+09 5.01e+12 Inf 1 0.034 0.9994
##
## Results are averaged over the levels of: Site
## P value adjustment: tukey method for comparing a family of 3 estimates
## Tests are performed on the log odds ratio scale- Check model assumptions
res <- resid(my_glmer1)
#produce residual vs. fitted plot
plot(fitted(my_glmer1), res)
#add a horizontal line at 0
abline(0,0)
#create Q-Q plot for residuals
qqnorm(res)
#add a straight diagonal line to the plot
qqline(res) 
2.3. Day 22
- Run LMM for survival at day 22
my_glmer1<- glmer(Survival_code ~ Site*Treatment + (1|Location), data=Juveniles_stat22_1, family = binomial)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## unable to evaluate scaled gradient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge: degenerate Hessian with 1 negative eigenvaluessummary(my_glmer1)## Warning in vcov.merMod(object, use.hessian = use.hessian): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX## Warning in vcov.merMod(object, correlation = correlation, sigm = sig): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Survival_code ~ Site * Treatment + (1 | Location)
## Data: Juveniles_stat22_1
##
## AIC BIC logLik deviance df.resid
## 74.6 96.2 -28.3 56.6 73
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.91312 0.00007 0.31820 0.43215 1.03905
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 0.8404 0.9167
## Number of obs: 82, groups: Location, 6
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.35408 0.93302 1.451 0.147
## SiteBonny Hills 17.31347 3498.86652 0.005 0.996
## SiteCronulla -0.06772 1.10986 -0.061 0.951
## SitePalm Beach 0.39728 1.20226 0.330 0.741
## TreatmentMHW1 0.44688 1.19831 0.373 0.709
## SiteBonny Hills:TreatmentMHW1 -18.02184 3498.86673 -0.005 0.996
## SiteCronulla:TreatmentMHW1 17.19573 3155.80978 0.005 0.996
## SitePalm Beach:TreatmentMHW1 0.42733 1.84602 0.231 0.817
##
## Correlation of Fixed Effects:
## (Intr) StBnnH StCrnl StPlmB TrMHW1 SBH:TM SC:TMH
## SitBnnyHlls 0.000
## SiteCronull -0.703 0.000
## SitePalmBch -0.652 0.000 0.550
## TretmntMHW1 -0.651 0.000 0.551 0.513
## StBHl:TMHW1 0.000 -1.000 0.000 0.000 0.000
## StCrn:TMHW1 0.000 0.000 0.000 0.000 0.000 0.000
## StPBc:TMHW1 0.428 0.000 -0.359 -0.651 -0.648 0.000 0.000
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## unable to evaluate scaled gradient
## Model failed to converge: degenerate Hessian with 1 negative eigenvaluesoverdisp.glmer(my_glmer1)## Residual deviance: 48.954 on 73 degrees of freedom (ratio: 0.671)my_glmer_anova <- car::Anova(my_glmer1)## Warning in vcov.merMod(mod, complete = FALSE): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RXmy_glmer_anova## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Survival_code
## Chisq Df Pr(>Chisq)
## Site 1.2830 3 0.7332
## Treatment 0.4715 1 0.4923
## Site:Treatment 0.0536 3 0.9967- Post-hoc test
emm3 <- emmeans(my_glmer1, specs = pairwise ~ Treatment, type = "response", adjust = "tukey")## Warning in vcov.merMod(object, complete = FALSE, ...): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX## NOTE: Results may be misleading due to involvement in interactionssummary(emm3)## $emmeans
## Treatment prob SE df asymp.LCL asymp.UCL
## Ambient 0.997 2.73 Inf 2.22e-16 1
## MHW1 0.998 1.74 Inf 2.22e-16 1
##
## Results are averaged over the levels of: Site
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scale
##
## $contrasts
## contrast odds.ratio SE df null z.ratio p.value
## Ambient / MHW1 0.707 832 Inf 1 0.000 0.9998
##
## Results are averaged over the levels of: Site
## Tests are performed on the log odds ratio scale- Check model assumptions
res <- resid(my_glmer1)
#produce residual vs. fitted plot
plot(fitted(my_glmer1), res)
#add a horizontal line at 0
abline(0,0)
#create Q-Q plot for residuals
qqnorm(res)
#add a straight diagonal line to the plot
qqline(res) 
3b) Structural intregrity loss
1. Visualise differences
Juv_means$Site <- factor(
Juv_means$Site,
levels = c("Bonny Hills","Black Head","Palm Beach","Cronulla"),
labels = c(
"Bonny Hills","Black Head","Palm Beach","Cronulla"
))
Juv_means$Day_equivalent1 <- as.numeric(Juv_means$Day_equivalent)
Melting_plot <- Juv_means %>%
filter(Day_equivalent == 8|Day_equivalent == 12|Day_equivalent == 18| Day_equivalent == 22) %>%
ggplot(aes(x= Day_equivalent, y= Mean_melting, colour = Site)) +
geom_errorbar(aes(ymin=Mean_melting-se_melting, ymax=Mean_melting+se_melting),colour = "grey", width=.2) +
geom_point() +
geom_line()+
Theme1 +
labs( y="Structural integrity loss (0-4)", x="Time (days)") +
scale_colour_manual(values = c("darkred", "orange", "darkolivegreen3", "royalblue3")) +
theme(strip.text.x = element_text(size = 12))+
scale_y_continuous(limits = c(0, 4.5), breaks = seq(0, 4, by = 1))+
scale_x_continuous(limits = c(7, 23), breaks = c(8, 12, 18,22))+
facet_wrap(~ factor(Treatment, c("Ambient","MHW1","MHW2")), ncol = 3)
Melting_plot 
ggsave("Melting_juv.jpeg", width = 220, height = 60, units = "mm")2. Statistical analyses
2.1. Day 8
- Run LMM for melting at day 8
my_lmer<- lmer(Melting ~ Site*Treatment + (1|Location), data=Juveniles_stat8_1)
summary(my_lmer)## Linear mixed model fit by REML ['lmerMod']
## Formula: Melting ~ Site * Treatment + (1 | Location)
## Data: Juveniles_stat8_1
##
## REML criterion at convergence: 201
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.9028 -0.6946 -0.2273 0.7249 3.3902
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 0.0623 0.2496
## Residual 0.2655 0.5153
## Number of obs: 121, groups: Location, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.55191 0.20095 2.747
## SiteBonny Hills -0.28746 0.24350 -1.181
## SiteCronulla -0.13525 0.22830 -0.592
## SitePalm Beach 0.02302 0.23313 0.099
## TreatmentMHW1 -0.27969 0.23861 -1.172
## TreatmentMHW2 0.37697 0.23813 1.583
## SiteBonny Hills:TreatmentMHW1 0.47657 0.34563 1.379
## SiteCronulla:TreatmentMHW1 0.11303 0.31810 0.355
## SitePalm Beach:TreatmentMHW1 0.15125 0.32539 0.465
## SiteBonny Hills:TreatmentMHW2 -0.30834 0.34482 -0.894
## SiteCronulla:TreatmentMHW2 -0.04364 0.31773 -0.137
## SitePalm Beach:TreatmentMHW2 -0.57634 0.33272 -1.732
##
## Correlation of Fixed Effects:
## (Intr) StBnnH StCrnl StPlmB TrMHW1 TrMHW2 SBH:TMHW1 SC:TMHW1
## SitBnnyHlls -0.609
## SiteCronull -0.654 0.536
## SitePalmBch -0.642 0.527 0.565
## TretmntMHW1 -0.628 0.513 0.553 0.542
## TretmntMHW2 -0.624 0.514 0.549 0.538 0.526
## StBHl:TMHW1 0.428 -0.700 -0.377 -0.368 -0.686 -0.359
## StCrn:TMHW1 0.471 -0.385 -0.719 -0.407 -0.750 -0.395 0.515
## StPBc:TMHW1 0.464 -0.379 -0.409 -0.720 -0.736 -0.389 0.503 0.552
## StBHl:TMHW2 0.428 -0.704 -0.377 -0.369 -0.360 -0.685 0.494 0.270
## StCrn:TMHW2 0.468 -0.385 -0.717 -0.403 -0.394 -0.749 0.269 0.515
## StPBc:TMHW2 0.451 -0.369 -0.397 -0.700 -0.381 -0.717 0.259 0.286
## SPB:TMHW1 SBH:TMHW2 SC:TMHW2
## SitBnnyHlls
## SiteCronull
## SitePalmBch
## TretmntMHW1
## TretmntMHW2
## StBHl:TMHW1
## StCrn:TMHW1
## StPBc:TMHW1
## StBHl:TMHW2 0.265
## StCrn:TMHW2 0.292 0.513
## StPBc:TMHW2 0.506 0.490 0.537my_glmer_anova <- car::Anova(my_lmer)
my_glmer_anova## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Melting
## Chisq Df Pr(>Chisq)
## Site 2.8005 3 0.4234
## Treatment 5.6570 2 0.0591 .
## Site:Treatment 9.6855 6 0.1385
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- Check model assumptions
res <- resid(my_lmer)
#produce residual vs. fitted plot
plot(fitted(my_lmer), res)
#add a horizontal line at 0
abline(0,0)
#create Q-Q plot for residuals
qqnorm(res)
#add a straight diagonal line to the plot
qqline(res) 
2.2. Day 12
- Run LMM for melting at day 12
my_glmer<- lmer(Melting ~ Site*Treatment + (1|Location), data=Juveniles_stat12_1)
summary(my_glmer)## Linear mixed model fit by REML ['lmerMod']
## Formula: Melting ~ Site * Treatment + (1 | Location)
## Data: Juveniles_stat12_1
##
## REML criterion at convergence: 354.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.6123 -0.6640 -0.1742 0.3761 3.0624
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 0.1835 0.4284
## Residual 1.0724 1.0355
## Number of obs: 122, groups: Location, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.26816 0.38938 3.257
## SiteBonny Hills -0.85786 0.48930 -1.753
## SiteCronulla -0.26816 0.45868 -0.585
## SitePalm Beach -0.52161 0.46833 -1.114
## TreatmentMHW1 -0.73542 0.46821 -1.571
## TreatmentMHW2 0.54142 0.47839 1.132
## SiteBonny Hills:TreatmentMHW1 1.45009 0.68672 2.112
## SiteCronulla:TreatmentMHW1 -0.01458 0.63083 -0.023
## SitePalm Beach:TreatmentMHW1 0.61912 0.64556 0.959
## SiteBonny Hills:TreatmentMHW2 -0.05573 0.69297 -0.080
## SiteCronulla:TreatmentMHW2 -0.12475 0.63842 -0.195
## SitePalm Beach:TreatmentMHW2 -0.21101 0.66849 -0.316
##
## Correlation of Fixed Effects:
## (Intr) StBnnH StCrnl StPlmB TrMHW1 TrMHW2 SBH:TMHW1 SC:TMHW1
## SitBnnyHlls -0.631
## SiteCronull -0.678 0.536
## SitePalmBch -0.665 0.527 0.565
## TretmntMHW1 -0.665 0.525 0.565 0.554
## TretmntMHW2 -0.647 0.514 0.549 0.538 0.538
## StBHl:TMHW1 0.448 -0.709 -0.380 -0.372 -0.678 -0.363
## StCrn:TMHW1 0.494 -0.390 -0.728 -0.411 -0.742 -0.400 0.503
## StPBc:TMHW1 0.486 -0.384 -0.413 -0.728 -0.728 -0.393 0.491 0.541
## StBHl:TMHW2 0.444 -0.704 -0.377 -0.369 -0.368 -0.685 0.500 0.273
## StCrn:TMHW2 0.485 -0.385 -0.717 -0.403 -0.403 -0.749 0.272 0.521
## StPBc:TMHW2 0.467 -0.369 -0.396 -0.700 -0.389 -0.717 0.261 0.289
## SPB:TMHW1 SBH:TMHW2 SC:TMHW2
## SitBnnyHlls
## SiteCronull
## SitePalmBch
## TretmntMHW1
## TretmntMHW2
## StBHl:TMHW1
## StCrn:TMHW1
## StPBc:TMHW1
## StBHl:TMHW2 0.268
## StCrn:TMHW2 0.295 0.513
## StPBc:TMHW2 0.511 0.491 0.537my_glmer_anova <- car::Anova(my_glmer)
my_glmer_anova## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Melting
## Chisq Df Pr(>Chisq)
## Site 2.6549 3 0.447940
## Treatment 10.4329 2 0.005427 **
## Site:Treatment 8.4187 6 0.209003
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1overdisp.glmer(my_glmer)## Residual deviance: 113.828 on 109 degrees of freedom (ratio: 1.044)- Post-hoc test
emm3 <- emmeans(my_glmer, specs = pairwise ~ Treatment, type = "response", adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionssummary(emm3)## $emmeans
## Treatment emmean SE df lower.CL upper.CL
## Ambient 0.856 0.240 10.3 0.324 1.39
## MHW1 0.634 0.238 10.2 0.104 1.16
## MHW2 1.300 0.244 10.9 0.762 1.84
##
## Results are averaged over the levels of: Site
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## Ambient - MHW1 0.222 0.230 105 0.963 0.6018
## Ambient - MHW2 -0.444 0.235 106 -1.886 0.1478
## MHW1 - MHW2 -0.665 0.234 106 -2.838 0.0149
##
## Results are averaged over the levels of: Site
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimates- Check model assumptions
res <- resid(my_glmer)
#produce residual vs. fitted plot
plot(fitted(my_glmer), res)
#add a horizontal line at 0
abline(0,0)
#create Q-Q plot for residuals
qqnorm(res)
#add a straight diagonal line to the plot
qqline(res) 
2.3. Day 18
- Run LMM for melting at day 18
my_glmer1<- lmer(Melting ~ Site*Treatment + (1|Location), data=Juveniles_stat18_1)
summary(my_glmer1)## Linear mixed model fit by REML ['lmerMod']
## Formula: Melting ~ Site * Treatment + (1 | Location)
## Data: Juveniles_stat18_1
##
## REML criterion at convergence: 313.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.80568 -0.55760 -0.04844 0.40106 2.76310
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 0.1388 0.3726
## Residual 1.1359 1.0658
## Number of obs: 107, groups: Location, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.264059 0.389187 3.248
## SiteBonny Hills -0.700787 0.503644 -1.391
## SiteCronulla -0.475431 0.490309 -0.970
## SitePalm Beach -0.247407 0.482187 -0.513
## TreatmentMHW1 -0.242245 0.493656 -0.491
## TreatmentMHW2 2.820516 0.564411 4.997
## SiteBonny Hills:TreatmentMHW1 0.845458 0.729836 1.158
## SiteCronulla:TreatmentMHW1 0.036950 0.670601 0.055
## SitePalm Beach:TreatmentMHW1 -0.002566 0.673776 -0.004
## SiteBonny Hills:TreatmentMHW2 0.116276 0.778468 0.149
## SiteCronulla:TreatmentMHW2 0.427020 0.757153 0.564
## SitePalm Beach:TreatmentMHW2 -0.341099 0.764936 -0.446
##
## Correlation of Fixed Effects:
## (Intr) StBnnH StCrnl StPlmB TrMHW1 TrMHW2 SBH:TMHW1 SC:TMHW1
## SitBnnyHlls -0.650
## SiteCronull -0.667 0.514
## SitePalmBch -0.685 0.527 0.539
## TretmntMHW1 -0.671 0.515 0.527 0.544
## TretmntMHW2 -0.580 0.448 0.460 0.469 0.459
## StBHl:TMHW1 0.449 -0.687 -0.353 -0.363 -0.672 -0.303
## StCrn:TMHW1 0.490 -0.377 -0.732 -0.397 -0.732 -0.337 0.492
## StPBc:TMHW1 0.495 -0.380 -0.387 -0.721 -0.737 -0.339 0.493 0.538
## StBHl:TMHW2 0.422 -0.647 -0.333 -0.341 -0.333 -0.722 0.444 0.244
## StCrn:TMHW2 0.433 -0.335 -0.648 -0.351 -0.344 -0.743 0.230 0.476
## StPBc:TMHW2 0.433 -0.332 -0.340 -0.631 -0.345 -0.735 0.230 0.251
## SPB:TMHW1 SBH:TMHW2 SC:TMHW2
## SitBnnyHlls
## SiteCronull
## SitePalmBch
## TretmntMHW1
## TretmntMHW2
## StBHl:TMHW1
## StCrn:TMHW1
## StPBc:TMHW1
## StBHl:TMHW2 0.246
## StCrn:TMHW2 0.254 0.538
## StPBc:TMHW2 0.456 0.533 0.549overdisp.glmer(my_glmer1)## Residual deviance: 104.119 on 94 degrees of freedom (ratio: 1.108)my_glmer_anova <- car::Anova(my_glmer1)
my_glmer_anova## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Melting
## Chisq Df Pr(>Chisq)
## Site 2.1805 3 0.5358
## Treatment 149.3260 2 <2e-16 ***
## Site:Treatment 3.6147 6 0.7286
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- Post-hoc test
emm3 <- emmeans(my_glmer1, specs = pairwise ~ Treatment, type = "response", adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionssummary(emm3)## $emmeans
## Treatment emmean SE df lower.CL upper.CL
## Ambient 0.908 0.231 11.4 0.403 1.41
## MHW1 0.886 0.231 11.5 0.381 1.39
## MHW2 3.779 0.260 16.5 3.229 4.33
##
## Results are averaged over the levels of: Site
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## Ambient - MHW1 0.0223 0.244 90.9 0.091 0.9954
## Ambient - MHW2 -2.8711 0.269 91.9 -10.675 <.0001
## MHW1 - MHW2 -2.8933 0.270 92.3 -10.701 <.0001
##
## Results are averaged over the levels of: Site
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimates- Check model assumptions
res <- resid(my_glmer1)
#produce residual vs. fitted plot
plot(fitted(my_glmer1), res)
#add a horizontal line at 0
abline(0,0)
#create Q-Q plot for residuals
qqnorm(res)
#add a straight diagonal line to the plot
qqline(res) 
2.4. Day 22
- Run LMM for melting at day 22
my_glmer1<- lmer(Melting ~ Site*Treatment + (1|Location), data=Juveniles_stat22_1)
summary(my_glmer1)## Linear mixed model fit by REML ['lmerMod']
## Formula: Melting ~ Site * Treatment + (1 | Location)
## Data: Juveniles_stat22_1
##
## REML criterion at convergence: 245.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8223 -0.7329 -0.1168 0.3833 2.2210
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 0.2738 0.5233
## Residual 1.4478 1.2032
## Number of obs: 77, groups: Location, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.46078 0.45863 3.185
## SiteBonny Hills -0.60545 0.56923 -1.064
## SiteCronulla -0.21399 0.55384 -0.386
## SitePalm Beach -0.25805 0.54587 -0.473
## TreatmentMHW1 -0.05691 0.57677 -0.099
## SiteBonny Hills:TreatmentMHW1 1.18794 0.83768 1.418
## SiteCronulla:TreatmentMHW1 -0.43988 0.77021 -0.571
## SitePalm Beach:TreatmentMHW1 0.29234 0.78338 0.373
##
## Correlation of Fixed Effects:
## (Intr) StBnnH StCrnl StPlmB TrMHW1 SBH:TM SC:TMH
## SitBnnyHlls -0.625
## SiteCronull -0.641 0.514
## SitePalmBch -0.660 0.528 0.539
## TretmntMHW1 -0.629 0.500 0.512 0.532
## StBHl:TMHW1 0.427 -0.676 -0.348 -0.359 -0.683
## StCrn:TMHW1 0.466 -0.372 -0.721 -0.393 -0.742 0.508
## StPBc:TMHW1 0.468 -0.373 -0.379 -0.705 -0.739 0.500 0.547overdisp.glmer(my_glmer1)## Residual deviance: 94.844 on 68 degrees of freedom (ratio: 1.395)my_glmer_anova <- car::Anova(my_glmer1)
my_glmer_anova## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Melting
## Chisq Df Pr(>Chisq)
## Site 1.7962 3 0.6158
## Treatment 0.2691 1 0.6040
## Site:Treatment 4.3079 3 0.2301- Post-hoc test
emm3 <- emmeans(my_glmer1, specs = pairwise ~ Treatment, type = "response", adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionssummary(emm3)## $emmeans
## Treatment emmean SE df lower.CL upper.CL
## Ambient 1.19 0.290 8.11 0.523 1.86
## MHW1 1.39 0.293 8.50 0.725 2.06
##
## Results are averaged over the levels of: Site
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## Ambient - MHW1 -0.203 0.281 65.1 -0.723 0.4722
##
## Results are averaged over the levels of: Site
## Degrees-of-freedom method: kenward-roger- Check model assumptions
res <- resid(my_glmer1)
#produce residual vs. fitted plot
plot(fitted(my_glmer1), res)
#add a horizontal line at 0
abline(0,0)
#create Q-Q plot for residuals
qqnorm(res)
#add a straight diagonal line to the plot
qqline(res) 
3c) Bleaching
1. Visualise differences
Bleaching_plot <- Juv_means %>%
filter(Day_equivalent == 8|Day_equivalent == 12|Day_equivalent == 18| Day_equivalent == 22) %>%
ggplot(aes(x= Day_equivalent, y= Mean_bleaching, colour = Site)) +
geom_errorbar(aes(ymin=Mean_bleaching-se_bleaching, ymax=Mean_bleaching+se_bleaching),colour = "grey", width=.2) +
geom_point() +
geom_line()+
Theme1 +
labs( y="Bleaching (0-4)", x="Time (days)") +
scale_colour_manual(values = c("darkred", "orange", "darkolivegreen3", "royalblue3")) +
theme(strip.text.x = element_text(size = 12))+
scale_y_continuous(limits = c(0, 4.5), breaks = seq(0, 4, by = 1))+
scale_x_continuous(limits = c(7, 23), breaks = c(8, 12, 18,22))+
facet_wrap(~ factor(Treatment, c("Ambient","MHW1","MHW2")), ncol = 3)
Bleaching_plot
ggsave("Bleaching_juv.jpeg", width = 220, height = 60, units = "mm")2. Statistical analyses
2.1. Day 8
- Run LMM for bleaching at day 8
Juveniles_stat8_1$log_Bleaching <- log(Juveniles_stat8_1$Bleaching +1)
my_lmer<- lmer(Bleaching ~ Site*Treatment + (1|Location), data=Juveniles_stat8_1)
summary(my_lmer)## Linear mixed model fit by REML ['lmerMod']
## Formula: Bleaching ~ Site * Treatment + (1 | Location)
## Data: Juveniles_stat8_1
##
## REML criterion at convergence: 126.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.62300 -0.52699 -0.27597 0.05107 2.51897
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 0.01147 0.1071
## Residual 0.13934 0.3733
## Number of obs: 121, groups: Location, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.43125 0.13262 3.252
## SiteBonny Hills -0.20698 0.17629 -1.174
## SiteCronulla -0.34792 0.16519 -2.106
## SitePalm Beach -0.23946 0.16860 -1.420
## TreatmentMHW1 -0.13340 0.17253 -0.773
## TreatmentMHW2 -0.34398 0.17225 -1.997
## SiteBonny Hills:TreatmentMHW1 0.05397 0.25021 0.216
## SiteCronulla:TreatmentMHW1 0.21674 0.23019 0.942
## SitePalm Beach:TreatmentMHW1 0.03002 0.23529 0.128
## SiteBonny Hills:TreatmentMHW2 0.75310 0.24976 3.015
## SiteCronulla:TreatmentMHW2 0.42731 0.22998 1.858
## SitePalm Beach:TreatmentMHW2 0.15424 0.24076 0.641
##
## Correlation of Fixed Effects:
## (Intr) StBnnH StCrnl StPlmB TrMHW1 TrMHW2 SBH:TMHW1 SC:TMHW1
## SitBnnyHlls -0.667
## SiteCronull -0.716 0.536
## SitePalmBch -0.702 0.526 0.564
## TretmntMHW1 -0.687 0.513 0.552 0.541
## TretmntMHW2 -0.684 0.514 0.549 0.538 0.526
## StBHl:TMHW1 0.469 -0.702 -0.377 -0.368 -0.687 -0.360
## StCrn:TMHW1 0.515 -0.384 -0.719 -0.405 -0.749 -0.394 0.515
## StPBc:TMHW1 0.507 -0.378 -0.407 -0.719 -0.736 -0.388 0.503 0.551
## StBHl:TMHW2 0.469 -0.704 -0.377 -0.369 -0.360 -0.685 0.495 0.270
## StCrn:TMHW2 0.512 -0.385 -0.717 -0.403 -0.394 -0.749 0.269 0.515
## StPBc:TMHW2 0.493 -0.369 -0.396 -0.700 -0.380 -0.716 0.259 0.285
## SPB:TMHW1 SBH:TMHW2 SC:TMHW2
## SitBnnyHlls
## SiteCronull
## SitePalmBch
## TretmntMHW1
## TretmntMHW2
## StBHl:TMHW1
## StCrn:TMHW1
## StPBc:TMHW1
## StBHl:TMHW2 0.265
## StCrn:TMHW2 0.291 0.513
## StPBc:TMHW2 0.505 0.491 0.537my_glmer_anova <- car::Anova(my_lmer)
my_glmer_anova## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Bleaching
## Chisq Df Pr(>Chisq)
## Site 7.0682 3 0.06975 .
## Treatment 0.2971 2 0.86194
## Site:Treatment 13.4601 6 0.03628 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- Post-hoc test
emm3 <- emmeans(my_lmer, specs = pairwise ~ Treatment|Site, type = "response", adjust = "tukey")
summary(emm3)## $emmeans
## Site = Black Head:
## Treatment emmean SE df lower.CL upper.CL
## Ambient 0.43125 0.133 76.2 0.1664 0.696
## MHW1 0.29785 0.126 71.6 0.0462 0.549
## MHW2 0.08727 0.127 70.4 -0.1658 0.340
##
## Site = Bonny Hills:
## Treatment emmean SE df lower.CL upper.CL
## Ambient 0.22427 0.132 77.2 -0.0393 0.488
## MHW1 0.14484 0.140 81.9 -0.1343 0.424
## MHW2 0.63339 0.140 81.8 0.3541 0.913
##
## Site = Cronulla:
## Treatment emmean SE df lower.CL upper.CL
## Ambient 0.08333 0.116 61.7 -0.1491 0.316
## MHW1 0.16667 0.116 61.7 -0.0658 0.399
## MHW2 0.16667 0.116 61.7 -0.0658 0.399
##
## Site = Palm Beach:
## Treatment emmean SE df lower.CL upper.CL
## Ambient 0.19179 0.121 66.4 -0.0497 0.433
## MHW1 0.08840 0.121 66.4 -0.1530 0.330
## MHW2 0.00205 0.132 77.2 -0.2615 0.266
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## Site = Black Head:
## contrast estimate SE df t.ratio p.value
## Ambient - MHW1 0.1334 0.173 105 0.771 0.7217
## Ambient - MHW2 0.3440 0.173 105 1.992 0.1191
## MHW1 - MHW2 0.2106 0.168 105 1.252 0.4257
##
## Site = Bonny Hills:
## contrast estimate SE df t.ratio p.value
## Ambient - MHW1 0.0794 0.182 105 0.436 0.9006
## Ambient - MHW2 -0.4091 0.182 105 -2.245 0.0684
## MHW1 - MHW2 -0.4885 0.188 105 -2.602 0.0284
##
## Site = Cronulla:
## contrast estimate SE df t.ratio p.value
## Ambient - MHW1 -0.0833 0.152 104 -0.547 0.8483
## Ambient - MHW2 -0.0833 0.152 104 -0.547 0.8483
## MHW1 - MHW2 0.0000 0.152 104 0.000 1.0000
##
## Site = Palm Beach:
## contrast estimate SE df t.ratio p.value
## Ambient - MHW1 0.1034 0.160 104 0.648 0.7939
## Ambient - MHW2 0.1897 0.168 104 1.129 0.4986
## MHW1 - MHW2 0.0864 0.168 104 0.514 0.8648
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimatesemm3 <- emmeans(my_lmer, specs = pairwise ~ Site|Treatment, type = "response", adjust = "tukey")
summary(emm3)## $emmeans
## Treatment = Ambient:
## Site emmean SE df lower.CL upper.CL
## Black Head 0.43125 0.133 76.2 0.1664 0.696
## Bonny Hills 0.22427 0.132 77.2 -0.0393 0.488
## Cronulla 0.08333 0.116 61.7 -0.1491 0.316
## Palm Beach 0.19179 0.121 66.4 -0.0497 0.433
##
## Treatment = MHW1:
## Site emmean SE df lower.CL upper.CL
## Black Head 0.29785 0.126 71.6 0.0462 0.549
## Bonny Hills 0.14484 0.140 81.9 -0.1343 0.424
## Cronulla 0.16667 0.116 61.7 -0.0658 0.399
## Palm Beach 0.08840 0.121 66.4 -0.1530 0.330
##
## Treatment = MHW2:
## Site emmean SE df lower.CL upper.CL
## Black Head 0.08727 0.127 70.4 -0.1658 0.340
## Bonny Hills 0.63339 0.140 81.8 0.3541 0.913
## Cronulla 0.16667 0.116 61.7 -0.0658 0.399
## Palm Beach 0.00205 0.132 77.2 -0.2615 0.266
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## Treatment = Ambient:
## contrast estimate SE df t.ratio p.value
## Black Head - Bonny Hills 0.2070 0.176 104 1.173 0.6452
## Black Head - Cronulla 0.3479 0.166 105 2.102 0.1592
## Black Head - Palm Beach 0.2395 0.169 105 1.417 0.4919
## Bonny Hills - Cronulla 0.1409 0.165 104 0.854 0.8282
## Bonny Hills - Palm Beach 0.0325 0.168 104 0.193 0.9974
## Cronulla - Palm Beach -0.1085 0.156 104 -0.695 0.8987
##
## Treatment = MHW1:
## contrast estimate SE df t.ratio p.value
## Black Head - Bonny Hills 0.1530 0.179 106 0.855 0.8278
## Black Head - Cronulla 0.1312 0.160 104 0.819 0.8452
## Black Head - Palm Beach 0.2094 0.164 105 1.280 0.5777
## Bonny Hills - Cronulla -0.0218 0.171 105 -0.127 0.9993
## Bonny Hills - Palm Beach 0.0564 0.175 105 0.323 0.9883
## Cronulla - Palm Beach 0.0783 0.156 104 0.502 0.9585
##
## Treatment = MHW2:
## contrast estimate SE df t.ratio p.value
## Black Head - Bonny Hills -0.5461 0.177 104 -3.078 0.0139
## Black Head - Cronulla -0.0794 0.161 105 -0.494 0.9602
## Black Head - Palm Beach 0.0852 0.172 105 0.495 0.9600
## Bonny Hills - Cronulla 0.4667 0.171 105 2.722 0.0375
## Bonny Hills - Palm Beach 0.6313 0.182 105 3.465 0.0042
## Cronulla - Palm Beach 0.1646 0.165 104 0.998 0.7510
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 4 estimates- Check model assumptions
res <- resid(my_lmer)
#produce residual vs. fitted plot
plot(fitted(my_lmer), res)
#add a horizontal line at 0
abline(0,0)
#create Q-Q plot for residuals
qqnorm(res)
#add a straight diagonal line to the plot
qqline(res) 
2.2. Day 12
- Run LMM for bleaching at day 12
my_glmer<- lmer(Bleaching ~ Site*Treatment + (1|Location), data=Juveniles_stat12_1)
summary(my_glmer)## Linear mixed model fit by REML ['lmerMod']
## Formula: Bleaching ~ Site * Treatment + (1 | Location)
## Data: Juveniles_stat12_1
##
## REML criterion at convergence: 123.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.7564 -0.5626 -0.2431 -0.0016 4.8662
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 0.001956 0.04423
## Residual 0.138680 0.37240
## Number of obs: 122, groups: Location, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.0005372 0.1257047 0.004
## SiteBonny Hills 0.1146815 0.1756600 0.653
## SiteCronulla 0.1661295 0.1644152 1.010
## SitePalm Beach 0.2750543 0.1676603 1.641
## TreatmentMHW1 0.0903600 0.1676571 0.539
## TreatmentMHW2 0.0006807 0.1713531 0.004
## SiteBonny Hills:TreatmentMHW1 0.0489268 0.2466286 0.198
## SiteCronulla:TreatmentMHW1 -0.0070266 0.2263234 -0.031
## SitePalm Beach:TreatmentMHW1 -0.1841903 0.2311108 -0.797
## SiteBonny Hills:TreatmentMHW2 0.1395448 0.2490874 0.560
## SiteCronulla:TreatmentMHW2 -0.0840140 0.2290748 -0.367
## SitePalm Beach:TreatmentMHW2 -0.2721646 0.2396454 -1.136
##
## Correlation of Fixed Effects:
## (Intr) StBnnH StCrnl StPlmB TrMHW1 TrMHW2 SBH:TMHW1 SC:TMHW1
## SitBnnyHlls -0.700
## SiteCronull -0.749 0.535
## SitePalmBch -0.735 0.525 0.562
## TretmntMHW1 -0.735 0.525 0.562 0.551
## TretmntMHW2 -0.718 0.513 0.549 0.538 0.538
## StBHl:TMHW1 0.498 -0.711 -0.381 -0.373 -0.679 -0.365
## StCrn:TMHW1 0.544 -0.389 -0.727 -0.408 -0.741 -0.399 0.503
## StPBc:TMHW1 0.534 -0.381 -0.408 -0.726 -0.726 -0.391 0.492 0.538
## StBHl:TMHW2 0.493 -0.705 -0.377 -0.370 -0.369 -0.686 0.501 0.274
## StCrn:TMHW2 0.537 -0.384 -0.717 -0.403 -0.403 -0.748 0.273 0.521
## StPBc:TMHW2 0.514 -0.367 -0.393 -0.700 -0.386 -0.715 0.261 0.286
## SPB:TMHW1 SBH:TMHW2 SC:TMHW2
## SitBnnyHlls
## SiteCronull
## SitePalmBch
## TretmntMHW1
## TretmntMHW2
## StBHl:TMHW1
## StCrn:TMHW1
## StPBc:TMHW1
## StBHl:TMHW2 0.268
## StCrn:TMHW2 0.293 0.513
## StPBc:TMHW2 0.508 0.491 0.535my_glmer_anova <- car::Anova(my_glmer)
my_glmer_anova## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Bleaching
## Chisq Df Pr(>Chisq)
## Site 3.5106 3 0.3194
## Treatment 1.8798 2 0.3907
## Site:Treatment 3.2022 6 0.7831overdisp.glmer(my_glmer)## Residual deviance: 15.103 on 109 degrees of freedom (ratio: 0.139)- Check model assumptions
res <- resid(my_glmer)
#produce residual vs. fitted plot
plot(fitted(my_glmer), res)
#add a horizontal line at 0
abline(0,0)
#create Q-Q plot for residuals
qqnorm(res)
#add a straight diagonal line to the plot
qqline(res) 
2.3. Day 18
- Run LMM for bleaching at day 18
my_glmer1<- lmer(Bleaching ~ Site*Treatment + (1|Location), data=Juveniles_stat18_1)
summary(my_glmer1)## Linear mixed model fit by REML ['lmerMod']
## Formula: Bleaching ~ Site * Treatment + (1 | Location)
## Data: Juveniles_stat18_1
##
## REML criterion at convergence: 177.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.3297 -0.7055 -0.3775 0.7989 2.7634
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 0.006089 0.07803
## Residual 0.284844 0.53371
## Number of obs: 107, groups: Location, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.33375 0.18128 1.841
## SiteBonny Hills 0.11475 0.25182 0.456
## SiteCronulla -0.13440 0.24534 -0.548
## SitePalm Beach 0.21562 0.24047 0.897
## TreatmentMHW1 0.17245 0.24596 0.701
## TreatmentMHW2 -0.17268 0.28179 -0.613
## SiteBonny Hills:TreatmentMHW1 -0.05235 0.36452 -0.144
## SiteCronulla:TreatmentMHW1 0.29487 0.33543 0.879
## SitePalm Beach:TreatmentMHW1 -0.44815 0.33561 -1.335
## SiteBonny Hills:TreatmentMHW2 0.15827 0.38943 0.406
## SiteCronulla:TreatmentMHW2 -0.03861 0.37871 -0.102
## SitePalm Beach:TreatmentMHW2 -0.09848 0.38222 -0.258
##
## Correlation of Fixed Effects:
## (Intr) StBnnH StCrnl StPlmB TrMHW1 TrMHW2 SBH:TMHW1 SC:TMHW1
## SitBnnyHlls -0.696
## SiteCronull -0.714 0.513
## SitePalmBch -0.731 0.525 0.538
## TretmntMHW1 -0.715 0.514 0.526 0.540
## TretmntMHW2 -0.622 0.448 0.459 0.469 0.459
## StBHl:TMHW1 0.481 -0.690 -0.354 -0.363 -0.673 -0.307
## StCrn:TMHW1 0.523 -0.376 -0.732 -0.395 -0.732 -0.336 0.493
## StPBc:TMHW1 0.526 -0.378 -0.386 -0.718 -0.735 -0.337 0.493 0.537
## StBHl:TMHW2 0.450 -0.647 -0.332 -0.340 -0.332 -0.723 0.446 0.243
## StCrn:TMHW2 0.463 -0.333 -0.648 -0.350 -0.342 -0.743 0.230 0.475
## StPBc:TMHW2 0.460 -0.330 -0.339 -0.629 -0.341 -0.736 0.229 0.249
## SPB:TMHW1 SBH:TMHW2 SC:TMHW2
## SitBnnyHlls
## SiteCronull
## SitePalmBch
## TretmntMHW1
## TretmntMHW2
## StBHl:TMHW1
## StCrn:TMHW1
## StPBc:TMHW1
## StBHl:TMHW2 0.244
## StCrn:TMHW2 0.252 0.537
## StPBc:TMHW2 0.452 0.533 0.548overdisp.glmer(my_glmer1)## Residual deviance: 26.682 on 94 degrees of freedom (ratio: 0.284)my_glmer_anova <- car::Anova(my_glmer1)
my_glmer_anova## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Bleaching
## Chisq Df Pr(>Chisq)
## Site 1.4252 3 0.69963
## Treatment 5.1810 2 0.07498 .
## Site:Treatment 7.0518 6 0.31608
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- Check model assumptions
res <- resid(my_glmer1)
#produce residual vs. fitted plot
plot(fitted(my_glmer1), res)
#add a horizontal line at 0
abline(0,0)
#create Q-Q plot for residuals
qqnorm(res)
#add a straight diagonal line to the plot
qqline(res) 
2.4. Day 22
- Run LMM for bleaching at day 22
my_glmer1<- lmer(Bleaching ~ Site*Treatment + (1|Location), data=Juveniles_stat22_1)
summary(my_glmer1)## Linear mixed model fit by REML ['lmerMod']
## Formula: Bleaching ~ Site * Treatment + (1 | Location)
## Data: Juveniles_stat22_1
##
## REML criterion at convergence: 140
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.6243 -0.7486 -0.4358 0.7251 2.6550
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 0.004843 0.06959
## Residual 0.338783 0.58205
## Number of obs: 77, groups: Location, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.447067 0.196546 2.275
## SiteBonny Hills 0.112237 0.274571 0.409
## SiteCronulla -0.146569 0.267529 -0.548
## SitePalm Beach 0.282819 0.262096 1.079
## TreatmentMHW1 0.002319 0.275320 0.008
## SiteBonny Hills:TreatmentMHW1 0.007431 0.402394 0.018
## SiteCronulla:TreatmentMHW1 0.613849 0.371062 1.654
## SitePalm Beach:TreatmentMHW1 -0.429557 0.375069 -1.145
##
## Correlation of Fixed Effects:
## (Intr) StBnnH StCrnl StPlmB TrMHW1 SBH:TM SC:TMH
## SitBnnyHlls -0.699
## SiteCronull -0.718 0.513
## SitePalmBch -0.735 0.525 0.538
## TretmntMHW1 -0.700 0.500 0.513 0.526
## StBHl:TMHW1 0.478 -0.682 -0.350 -0.358 -0.683
## StCrn:TMHW1 0.518 -0.371 -0.721 -0.389 -0.741 0.506
## StPBc:TMHW1 0.515 -0.368 -0.377 -0.700 -0.734 0.501 0.544overdisp.glmer(my_glmer1)## Residual deviance: 23.121 on 68 degrees of freedom (ratio: 0.34)my_glmer_anova <- car::Anova(my_glmer1)
my_glmer_anova## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Bleaching
## Chisq Df Pr(>Chisq)
## Site 1.0441 3 0.79058
## Treatment 0.2085 1 0.64792
## Site:Treatment 8.7135 3 0.03335 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- Post-hoc test
emm3 <- emmeans(my_glmer1, specs = pairwise ~ Treatment|Site, type = "response", adjust = "tukey")
summary(emm3)## $emmeans
## Site = Black Head:
## Treatment emmean SE df lower.CL upper.CL
## Ambient 0.447 0.199 61.2 0.0501 0.844
## MHW1 0.449 0.199 61.4 0.0522 0.847
##
## Site = Bonny Hills:
## Treatment emmean SE df lower.CL upper.CL
## Ambient 0.559 0.197 63.4 0.1653 0.953
## MHW1 0.569 0.227 62.7 0.1158 1.022
##
## Site = Cronulla:
## Treatment emmean SE df lower.CL upper.CL
## Ambient 0.300 0.187 61.5 -0.0735 0.674
## MHW1 0.917 0.170 58.3 0.5756 1.258
##
## Site = Palm Beach:
## Treatment emmean SE df lower.CL upper.CL
## Ambient 0.730 0.178 59.9 0.3734 1.086
## MHW1 0.303 0.187 61.6 -0.0713 0.677
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## Site = Black Head:
## contrast estimate SE df t.ratio p.value
## Ambient - MHW1 -0.00232 0.279 67.9 -0.008 0.9934
##
## Site = Bonny Hills:
## contrast estimate SE df t.ratio p.value
## Ambient - MHW1 -0.00975 0.296 66.3 -0.033 0.9738
##
## Site = Cronulla:
## contrast estimate SE df t.ratio p.value
## Ambient - MHW1 -0.61617 0.250 64.7 -2.466 0.0163
##
## Site = Palm Beach:
## contrast estimate SE df t.ratio p.value
## Ambient - MHW1 0.42724 0.255 65.3 1.672 0.0992
##
## Degrees-of-freedom method: kenward-rogeremm3 <- emmeans(my_glmer1, specs = pairwise ~ Site|Treatment, type = "response", adjust = "tukey")
summary(emm3)## $emmeans
## Treatment = Ambient:
## Site emmean SE df lower.CL upper.CL
## Black Head 0.447 0.199 61.2 0.0501 0.844
## Bonny Hills 0.559 0.197 63.4 0.1653 0.953
## Cronulla 0.300 0.187 61.5 -0.0735 0.674
## Palm Beach 0.730 0.178 59.9 0.3734 1.086
##
## Treatment = MHW1:
## Site emmean SE df lower.CL upper.CL
## Black Head 0.449 0.199 61.4 0.0522 0.847
## Bonny Hills 0.569 0.227 62.7 0.1158 1.022
## Cronulla 0.917 0.170 58.3 0.5756 1.258
## Palm Beach 0.303 0.187 61.6 -0.0713 0.677
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## Treatment = Ambient:
## contrast estimate SE df t.ratio p.value
## Black Head - Bonny Hills -0.112 0.275 65.1 -0.408 0.9769
## Black Head - Cronulla 0.147 0.268 64.6 0.547 0.9470
## Black Head - Palm Beach -0.283 0.264 66.5 -1.071 0.7085
## Bonny Hills - Cronulla 0.259 0.268 64.5 0.966 0.7690
## Bonny Hills - Palm Beach -0.171 0.262 64.7 -0.650 0.9150
## Cronulla - Palm Beach -0.429 0.255 65.3 -1.681 0.3420
##
## Treatment = MHW1:
## contrast estimate SE df t.ratio p.value
## Black Head - Bonny Hills -0.120 0.299 68.2 -0.400 0.9782
## Black Head - Cronulla -0.467 0.259 65.8 -1.807 0.2795
## Black Head - Palm Beach 0.147 0.269 65.5 0.546 0.9474
## Bonny Hills - Cronulla -0.348 0.281 67.3 -1.238 0.6053
## Bonny Hills - Palm Beach 0.266 0.293 68.3 0.909 0.8000
## Cronulla - Palm Beach 0.614 0.250 64.7 2.458 0.0765
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 4 estimates- Check model assumptions
res <- resid(my_glmer1)
#produce residual vs. fitted plot
plot(fitted(my_glmer1), res)
#add a horizontal line at 0
abline(0,0)
#create Q-Q plot for residuals
qqnorm(res)
#add a straight diagonal line to the plot
qqline(res) 
3d) Fouling
1. Visualise differences
Fouling_plot <- Juv_means %>%
filter(Day_equivalent == 8|Day_equivalent == 12|Day_equivalent == 18| Day_equivalent == 22) %>%
ggplot(aes(x= Day_equivalent, y= Mean_fouling, colour = Site)) +
geom_errorbar(aes(ymin=Mean_fouling-se_fouling, ymax=Mean_fouling+se_fouling),colour = "grey", width=.2) +
geom_point() +
geom_line()+
Theme1 +
labs( y="Fouling (0-4)", x="Time (days)") +
scale_colour_manual(values = c("darkred", "orange", "darkolivegreen3", "royalblue3")) +
theme(strip.text.x = element_text(size = 12))+
scale_y_continuous(limits = c(0, 4.5), breaks = seq(0, 4, by = 1))+
scale_x_continuous(limits = c(7, 23), breaks = c(8, 12, 18,22))+
facet_wrap(~ factor(Treatment, c("Ambient","MHW1","MHW2")), ncol = 3)
Fouling_plot
ggsave("Fouling_juv.jpeg", width = 220, height = 60, units = "mm")2. Statistical analyses
2.1. Day 8
- Run LMM for fouling at day 8
my_lmer<- lmer(Fouling ~ Site*Treatment + (1|Location), data=Juveniles_stat8_1)## boundary (singular) fit: see help('isSingular')summary(my_lmer)## Linear mixed model fit by REML ['lmerMod']
## Formula: Fouling ~ Site * Treatment + (1 | Location)
## Data: Juveniles_stat8_1
##
## REML criterion at convergence: 244.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.5254 -0.9707 0.3051 0.5547 2.0800
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 0.0000 0.0000
## Residual 0.4298 0.6556
## Number of obs: 121, groups: Location, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.000e+00 2.185e-01 4.576
## SiteBonny Hills -1.047e-15 3.090e-01 0.000
## SiteCronulla -3.333e-01 2.891e-01 -1.153
## SitePalm Beach -3.636e-01 2.947e-01 -1.234
## TreatmentMHW1 -2.000e-01 3.012e-01 -0.664
## TreatmentMHW2 -3.000e-01 3.012e-01 -0.996
## SiteBonny Hills:TreatmentMHW1 -3.000e-01 4.384e-01 -0.684
## SiteCronulla:TreatmentMHW1 3.333e-02 4.029e-01 0.083
## SitePalm Beach:TreatmentMHW1 2.000e-01 4.109e-01 0.487
## SiteBonny Hills:TreatmentMHW2 5.000e-02 4.384e-01 0.114
## SiteCronulla:TreatmentMHW2 3.000e-01 4.029e-01 0.745
## SitePalm Beach:TreatmentMHW2 2.192e-01 4.214e-01 0.520
##
## Correlation of Fixed Effects:
## (Intr) StBnnH StCrnl StPlmB TrMHW1 TrMHW2 SBH:TMHW1 SC:TMHW1
## SitBnnyHlls -0.707
## SiteCronull -0.756 0.535
## SitePalmBch -0.742 0.524 0.561
## TretmntMHW1 -0.725 0.513 0.548 0.538
## TretmntMHW2 -0.725 0.513 0.548 0.538 0.526
## StBHl:TMHW1 0.498 -0.705 -0.377 -0.370 -0.687 -0.362
## StCrn:TMHW1 0.542 -0.383 -0.717 -0.402 -0.748 -0.393 0.514
## StPBc:TMHW1 0.532 -0.376 -0.402 -0.717 -0.733 -0.386 0.504 0.548
## StBHl:TMHW2 0.498 -0.705 -0.377 -0.370 -0.362 -0.687 0.497 0.270
## StCrn:TMHW2 0.542 -0.383 -0.717 -0.402 -0.393 -0.748 0.270 0.515
## StPBc:TMHW2 0.519 -0.367 -0.392 -0.699 -0.376 -0.715 0.258 0.281
## SPB:TMHW1 SBH:TMHW2 SC:TMHW2
## SitBnnyHlls
## SiteCronull
## SitePalmBch
## TretmntMHW1
## TretmntMHW2
## StBHl:TMHW1
## StCrn:TMHW1
## StPBc:TMHW1
## StBHl:TMHW2 0.265
## StCrn:TMHW2 0.288 0.514
## StPBc:TMHW2 0.501 0.491 0.534
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')my_glmer_anova <- car::Anova(my_lmer)
my_glmer_anova## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Fouling
## Chisq Df Pr(>Chisq)
## Site 2.5669 3 0.4633
## Treatment 1.9772 2 0.3721
## Site:Treatment 2.1223 6 0.9081- Check model assumptions
res <- resid(my_lmer)
#produce residual vs. fitted plot
plot(fitted(my_lmer), res)
#add a horizontal line at 0
abline(0,0)
#create Q-Q plot for residuals
qqnorm(res)
#add a straight diagonal line to the plot
qqline(res) 
2.2. Day 12
- Run LMM for fouling at day 12
my_glmer<- lmer(Fouling ~ Site*Treatment + (1|Location), data=Juveniles_stat12_1)## boundary (singular) fit: see help('isSingular')summary(my_glmer)## Linear mixed model fit by REML ['lmerMod']
## Formula: Fouling ~ Site * Treatment + (1 | Location)
## Data: Juveniles_stat12_1
##
## REML criterion at convergence: 208.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.6172 -0.8086 -0.3032 0.8270 2.6464
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 0.0000 0.0000
## Residual 0.3021 0.5496
## Number of obs: 122, groups: Location, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.66667 0.18321 3.639
## SiteBonny Hills 0.22222 0.25910 0.858
## SiteCronulla -0.16667 0.24237 -0.688
## SitePalm Beach -0.12121 0.24704 -0.491
## TreatmentMHW1 0.06061 0.24704 0.245
## TreatmentMHW2 -0.26667 0.25254 -1.056
## SiteBonny Hills:TreatmentMHW1 -0.57449 0.36381 -1.579
## SiteCronulla:TreatmentMHW1 -0.39394 0.33374 -1.180
## SitePalm Beach:TreatmentMHW1 -0.06061 0.34052 -0.178
## SiteBonny Hills:TreatmentMHW2 -0.24722 0.36757 -0.673
## SiteCronulla:TreatmentMHW2 -0.06667 0.33783 -0.197
## SitePalm Beach:TreatmentMHW2 0.16566 0.35328 0.469
##
## Correlation of Fixed Effects:
## (Intr) StBnnH StCrnl StPlmB TrMHW1 TrMHW2 SBH:TMHW1 SC:TMHW1
## SitBnnyHlls -0.707
## SiteCronull -0.756 0.535
## SitePalmBch -0.742 0.524 0.561
## TretmntMHW1 -0.742 0.524 0.561 0.550
## TretmntMHW2 -0.725 0.513 0.548 0.538 0.538
## StBHl:TMHW1 0.504 -0.712 -0.381 -0.373 -0.679 -0.365
## StCrn:TMHW1 0.549 -0.388 -0.726 -0.407 -0.740 -0.398 0.503
## StPBc:TMHW1 0.538 -0.380 -0.407 -0.725 -0.725 -0.390 0.493 0.537
## StBHl:TMHW2 0.498 -0.705 -0.377 -0.370 -0.370 -0.687 0.502 0.274
## StCrn:TMHW2 0.542 -0.383 -0.717 -0.402 -0.402 -0.748 0.273 0.521
## StPBc:TMHW2 0.519 -0.367 -0.392 -0.699 -0.385 -0.715 0.261 0.285
## SPB:TMHW1 SBH:TMHW2 SC:TMHW2
## SitBnnyHlls
## SiteCronull
## SitePalmBch
## TretmntMHW1
## TretmntMHW2
## StBHl:TMHW1
## StCrn:TMHW1
## StPBc:TMHW1
## StBHl:TMHW2 0.268
## StCrn:TMHW2 0.292 0.514
## StPBc:TMHW2 0.507 0.491 0.534
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')my_glmer_anova <- car::Anova(my_glmer)
my_glmer_anova## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Fouling
## Chisq Df Pr(>Chisq)
## Site 6.9582 3 0.07324 .
## Treatment 6.0994 2 0.04737 *
## Site:Treatment 4.1402 6 0.65771
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1overdisp.glmer(my_glmer)## Residual deviance: 33.231 on 109 degrees of freedom (ratio: 0.305)- Post-hoc test
emm3 <- emmeans(my_glmer, specs = pairwise ~ Treatment, type = "response", adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionssummary(emm3)## $emmeans
## Treatment emmean SE df lower.CL upper.CL
## Ambient 0.650 0.0870 35.4 0.474 0.827
## MHW1 0.454 0.0860 35.9 0.279 0.628
## MHW2 0.347 0.0898 36.2 0.164 0.529
##
## Results are averaged over the levels of: Site
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## Ambient - MHW1 0.197 0.122 106 1.610 0.2460
## Ambient - MHW2 0.304 0.125 107 2.438 0.0430
## MHW1 - MHW2 0.107 0.124 107 0.863 0.6650
##
## Results are averaged over the levels of: Site
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimates- Check model assumptions
res <- resid(my_glmer)
#produce residual vs. fitted plot
plot(fitted(my_glmer), res)
#add a horizontal line at 0
abline(0,0)
#create Q-Q plot for residuals
qqnorm(res)
#add a straight diagonal line to the plot
qqline(res) 
2.3 Day 18
- Run LMM for fouling at day 18
my_glmer1<- lmer(Fouling ~ Site*Treatment + (1|Location), data=Juveniles_stat18_1)## boundary (singular) fit: see help('isSingular')summary(my_glmer1)## Linear mixed model fit by REML ['lmerMod']
## Formula: Fouling ~ Site * Treatment + (1 | Location)
## Data: Juveniles_stat18_1
##
## REML criterion at convergence: 142
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.4960 -0.4987 -0.1870 0.0000 3.4679
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 0.0000 0.0000
## Residual 0.1986 0.4456
## Number of obs: 107, groups: Location, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.22222 0.14855 1.496
## SiteBonny Hills 0.44444 0.21008 2.116
## SiteCronulla -0.02222 0.20476 -0.109
## SitePalm Beach 0.23232 0.20030 1.160
## TreatmentMHW1 -0.02222 0.20476 -0.109
## TreatmentMHW2 -0.22222 0.23488 -0.946
## SiteBonny Hills:TreatmentMHW1 -0.35873 0.30392 -1.180
## SiteCronulla:TreatmentMHW1 -0.09444 0.27989 -0.337
## SitePalm Beach:TreatmentMHW1 -0.15960 0.27935 -0.571
## SiteBonny Hills:TreatmentMHW2 -0.44444 0.32497 -1.368
## SiteCronulla:TreatmentMHW2 0.02222 0.31599 0.070
## SitePalm Beach:TreatmentMHW2 -0.08947 0.31874 -0.281
##
## Correlation of Fixed Effects:
## (Intr) StBnnH StCrnl StPlmB TrMHW1 TrMHW2 SBH:TMHW1 SC:TMHW1
## SitBnnyHlls -0.707
## SiteCronull -0.725 0.513
## SitePalmBch -0.742 0.524 0.538
## TretmntMHW1 -0.725 0.513 0.526 0.538
## TretmntMHW2 -0.632 0.447 0.459 0.469 0.459
## StBHl:TMHW1 0.489 -0.691 -0.355 -0.362 -0.674 -0.309
## StCrn:TMHW1 0.531 -0.375 -0.732 -0.394 -0.732 -0.336 0.493
## StPBc:TMHW1 0.532 -0.376 -0.386 -0.717 -0.733 -0.336 0.494 0.536
## StBHl:TMHW2 0.457 -0.646 -0.332 -0.339 -0.332 -0.723 0.447 0.243
## StCrn:TMHW2 0.470 -0.332 -0.648 -0.349 -0.341 -0.743 0.230 0.474
## StPBc:TMHW2 0.466 -0.330 -0.338 -0.628 -0.338 -0.737 0.228 0.247
## SPB:TMHW1 SBH:TMHW2 SC:TMHW2
## SitBnnyHlls
## SiteCronull
## SitePalmBch
## TretmntMHW1
## TretmntMHW2
## StBHl:TMHW1
## StCrn:TMHW1
## StPBc:TMHW1
## StBHl:TMHW2 0.243
## StCrn:TMHW2 0.250 0.537
## StPBc:TMHW2 0.451 0.533 0.548
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')overdisp.glmer(my_glmer1)## Residual deviance: 18.867 on 94 degrees of freedom (ratio: 0.201)my_glmer_anova <- car::Anova(my_glmer1)
my_glmer_anova## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Fouling
## Chisq Df Pr(>Chisq)
## Site 5.4415 3 0.142179
## Treatment 9.9190 2 0.007016 **
## Site:Treatment 3.2595 6 0.775641
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- Post-hoc test
emm3 <- emmeans(my_glmer1, specs = pairwise ~ Treatment, type = "response", adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionssummary(emm3)## $emmeans
## Treatment emmean SE df lower.CL upper.CL
## Ambient 0.3859 0.0723 28.9 0.2380 0.534
## MHW1 0.2104 0.0723 30.6 0.0628 0.358
## MHW2 0.0357 0.0866 36.8 -0.1399 0.211
##
## Results are averaged over the levels of: Site
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## Ambient - MHW1 0.175 0.102 92.2 1.720 0.2031
## Ambient - MHW2 0.350 0.112 93.6 3.134 0.0065
## MHW1 - MHW2 0.175 0.112 94.1 1.559 0.2687
##
## Results are averaged over the levels of: Site
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimates- Check model assumptions
res <- resid(my_glmer1)
#produce residual vs. fitted plot
plot(fitted(my_glmer1), res)
#add a horizontal line at 0
abline(0,0)
#create Q-Q plot for residuals
qqnorm(res)
#add a straight diagonal line to the plot
qqline(res) 
2.3. Day 22
- Run LMM for fouling at day 22
my_glmer1<- lmer(Fouling ~ Site*Treatment + (1|Location), data=Juveniles_stat22_1)## boundary (singular) fit: see help('isSingular')summary(my_glmer1)## Linear mixed model fit by REML ['lmerMod']
## Formula: Fouling ~ Site * Treatment + (1 | Location)
## Data: Juveniles_stat22_1
##
## REML criterion at convergence: 110.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.9587 -0.4687 -0.2344 0.0000 3.2809
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 0.0000 0.0000
## Residual 0.2248 0.4741
## Number of obs: 77, groups: Location, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.11111 0.15804 0.703
## SiteBonny Hills 0.33333 0.22350 1.491
## SiteCronulla 0.08889 0.21784 0.408
## SitePalm Beach 0.34343 0.21310 1.612
## TreatmentMHW1 0.11111 0.22350 0.497
## SiteBonny Hills:TreatmentMHW1 -0.55556 0.32718 -1.698
## SiteCronulla:TreatmentMHW1 -0.22778 0.30194 -0.754
## SitePalm Beach:TreatmentMHW1 -0.36566 0.30474 -1.200
##
## Correlation of Fixed Effects:
## (Intr) StBnnH StCrnl StPlmB TrMHW1 SBH:TM SC:TMH
## SitBnnyHlls -0.707
## SiteCronull -0.725 0.513
## SitePalmBch -0.742 0.524 0.538
## TretmntMHW1 -0.707 0.500 0.513 0.524
## StBHl:TMHW1 0.483 -0.683 -0.350 -0.358 -0.683
## StCrn:TMHW1 0.523 -0.370 -0.721 -0.388 -0.740 0.506
## StPBc:TMHW1 0.519 -0.367 -0.376 -0.699 -0.733 0.501 0.543
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')overdisp.glmer(my_glmer1)## Residual deviance: 15.511 on 68 degrees of freedom (ratio: 0.228)my_glmer_anova <- car::Anova(my_glmer1)
my_glmer_anova## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Fouling
## Chisq Df Pr(>Chisq)
## Site 1.9312 3 0.5868
## Treatment 2.4106 1 0.1205
## Site:Treatment 3.1366 3 0.3710- Post-hoc test
emm3 <- emmeans(my_glmer1, specs = pairwise ~ Treatment|Site, type = "response", adjust = "tukey")
summary(emm3)## $emmeans
## Site = Black Head:
## Treatment emmean SE df lower.CL upper.CL
## Ambient 0.1111 0.160 63.2 -0.2084 0.431
## MHW1 0.2222 0.160 63.2 -0.0977 0.542
##
## Site = Bonny Hills:
## Treatment emmean SE df lower.CL upper.CL
## Ambient 0.4444 0.159 65.4 0.1271 0.762
## MHW1 0.0000 0.183 63.6 -0.3655 0.365
##
## Site = Cronulla:
## Treatment emmean SE df lower.CL upper.CL
## Ambient 0.2000 0.151 64.2 -0.1007 0.501
## MHW1 0.0833 0.137 62.0 -0.1903 0.357
##
## Site = Palm Beach:
## Treatment emmean SE df lower.CL upper.CL
## Ambient 0.4545 0.143 63.0 0.1682 0.741
## MHW1 0.2000 0.151 64.2 -0.1008 0.501
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## Site = Black Head:
## contrast estimate SE df t.ratio p.value
## Ambient - MHW1 -0.111 0.227 68.2 -0.489 0.6264
##
## Site = Bonny Hills:
## contrast estimate SE df t.ratio p.value
## Ambient - MHW1 0.444 0.241 66.5 1.845 0.0695
##
## Site = Cronulla:
## contrast estimate SE df t.ratio p.value
## Ambient - MHW1 0.117 0.203 64.8 0.573 0.5683
##
## Site = Palm Beach:
## contrast estimate SE df t.ratio p.value
## Ambient - MHW1 0.255 0.208 65.4 1.224 0.2254
##
## Degrees-of-freedom method: kenward-rogeremm3 <- emmeans(my_glmer1, specs = pairwise ~ Site|Treatment, type = "response", adjust = "tukey")
summary(emm3)## $emmeans
## Treatment = Ambient:
## Site emmean SE df lower.CL upper.CL
## Black Head 0.1111 0.160 63.2 -0.2084 0.431
## Bonny Hills 0.4444 0.159 65.4 0.1271 0.762
## Cronulla 0.2000 0.151 64.2 -0.1007 0.501
## Palm Beach 0.4545 0.143 63.0 0.1682 0.741
##
## Treatment = MHW1:
## Site emmean SE df lower.CL upper.CL
## Black Head 0.2222 0.160 63.2 -0.0977 0.542
## Bonny Hills 0.0000 0.183 63.6 -0.3655 0.365
## Cronulla 0.0833 0.137 62.0 -0.1903 0.357
## Palm Beach 0.2000 0.151 64.2 -0.1008 0.501
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## Treatment = Ambient:
## contrast estimate SE df t.ratio p.value
## Black Head - Bonny Hills -0.3333 0.224 65.2 -1.487 0.4512
## Black Head - Cronulla -0.0889 0.218 64.7 -0.407 0.9770
## Black Head - Palm Beach -0.3434 0.215 66.7 -1.598 0.3868
## Bonny Hills - Cronulla 0.2444 0.218 64.6 1.121 0.6783
## Bonny Hills - Palm Beach -0.0101 0.214 64.9 -0.047 1.0000
## Cronulla - Palm Beach -0.2545 0.208 65.4 -1.224 0.6141
##
## Treatment = MHW1:
## contrast estimate SE df t.ratio p.value
## Black Head - Bonny Hills 0.2222 0.244 68.5 0.913 0.7982
## Black Head - Cronulla 0.1389 0.211 66.1 0.659 0.9119
## Black Head - Palm Beach 0.0222 0.219 65.7 0.102 0.9996
## Bonny Hills - Cronulla -0.0833 0.228 67.6 -0.365 0.9833
## Bonny Hills - Palm Beach -0.2000 0.238 68.6 -0.840 0.8353
## Cronulla - Palm Beach -0.1167 0.203 64.8 -0.573 0.9397
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 4 estimates- Check model assumptions
res <- resid(my_glmer1)
#produce residual vs. fitted plot
plot(fitted(my_glmer1), res)
#add a horizontal line at 0
abline(0,0)
#create Q-Q plot for residuals
qqnorm(res)
#add a straight diagonal line to the plot
qqline(res) 
3e) Green colouration
1. Visualise differences
Green_plot <- Juv_means %>%
filter(Day_equivalent == 8|Day_equivalent == 12|Day_equivalent == 18| Day_equivalent == 22) %>%
ggplot(aes(x= Day_equivalent, y= Mean_green, colour = Site)) +
geom_errorbar(aes(ymin=Mean_green-se_green, ymax=Mean_green+se_green),colour = "grey", width=.2) +
geom_point() +
geom_line()+
Theme1 +
labs( y="Green colouration (0-4)", x="Time (days)") +
scale_colour_manual(values = c("darkred", "orange", "darkolivegreen3", "royalblue3")) +
theme(strip.text.x = element_text(size = 12))+
scale_y_continuous(limits = c(0, 4.5), breaks = seq(0, 4, by = 1))+
scale_x_continuous(limits = c(7, 23), breaks = c(8, 12, 18,22))+
facet_wrap(~ factor(Treatment, c("Ambient","MHW1","MHW2")), ncol = 3)
Green_plot
ggsave("Green_juv.jpeg", width = 220, height = 60, units = "mm")2. Statistical analyses
2.1. Day 8
- Run LMM for green colouration at day 8
my_lmer<- lmer(Green ~ Site*Treatment + (1|Location), data=Juveniles_stat8_1)## boundary (singular) fit: see help('isSingular')summary(my_lmer)## Linear mixed model fit by REML ['lmerMod']
## Formula: Green ~ Site * Treatment + (1 | Location)
## Data: Juveniles_stat8_1
##
## REML criterion at convergence: 29.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.047 0.000 0.000 0.000 6.143
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 0.00000 0.0000
## Residual 0.05963 0.2442
## Number of obs: 121, groups: Location, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -1.238e-16 8.140e-02 0.000
## SiteBonny Hills 1.773e-16 1.151e-01 0.000
## SiteCronulla 1.404e-16 1.077e-01 0.000
## SitePalm Beach 0.000e+00 1.098e-01 0.000
## TreatmentMHW1 1.173e-16 1.122e-01 0.000
## TreatmentMHW2 5.000e-01 1.122e-01 4.456
## SiteBonny Hills:TreatmentMHW1 -1.598e-16 1.633e-01 0.000
## SiteCronulla:TreatmentMHW1 -9.887e-17 1.501e-01 0.000
## SitePalm Beach:TreatmentMHW1 6.018e-17 1.531e-01 0.000
## SiteBonny Hills:TreatmentMHW2 -5.000e-01 1.633e-01 -3.062
## SiteCronulla:TreatmentMHW2 -5.000e-01 1.501e-01 -3.331
## SitePalm Beach:TreatmentMHW2 -5.000e-01 1.570e-01 -3.186
##
## Correlation of Fixed Effects:
## (Intr) StBnnH StCrnl StPlmB TrMHW1 TrMHW2 SBH:TMHW1 SC:TMHW1
## SitBnnyHlls -0.707
## SiteCronull -0.756 0.535
## SitePalmBch -0.742 0.524 0.561
## TretmntMHW1 -0.725 0.513 0.548 0.538
## TretmntMHW2 -0.725 0.513 0.548 0.538 0.526
## StBHl:TMHW1 0.498 -0.705 -0.377 -0.370 -0.687 -0.362
## StCrn:TMHW1 0.542 -0.383 -0.717 -0.402 -0.748 -0.393 0.514
## StPBc:TMHW1 0.532 -0.376 -0.402 -0.717 -0.733 -0.386 0.504 0.548
## StBHl:TMHW2 0.498 -0.705 -0.377 -0.370 -0.362 -0.687 0.497 0.270
## StCrn:TMHW2 0.542 -0.383 -0.717 -0.402 -0.393 -0.748 0.270 0.515
## StPBc:TMHW2 0.519 -0.367 -0.392 -0.699 -0.376 -0.715 0.258 0.281
## SPB:TMHW1 SBH:TMHW2 SC:TMHW2
## SitBnnyHlls
## SiteCronull
## SitePalmBch
## TretmntMHW1
## TretmntMHW2
## StBHl:TMHW1
## StCrn:TMHW1
## StPBc:TMHW1
## StBHl:TMHW2 0.265
## StCrn:TMHW2 0.288 0.514
## StPBc:TMHW2 0.501 0.491 0.534
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')my_glmer_anova <- car::Anova(my_lmer)
my_glmer_anova## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Green
## Chisq Df Pr(>Chisq)
## Site 10.5346 3 0.01453 *
## Treatment 6.8278 2 0.03291 *
## Site:Treatment 20.6390 6 0.00213 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- Post-hoc test
emm3 <- emmeans(my_lmer, specs = pairwise ~ Treatment|Site, type = "response", adjust = "tukey")
summary(emm3)## $emmeans
## Site = Black Head:
## Treatment emmean SE df lower.CL upper.CL
## Ambient 0.0 0.0820 103 -0.163 0.163
## MHW1 0.0 0.0774 104 -0.154 0.154
## MHW2 0.5 0.0777 102 0.346 0.654
##
## Site = Bonny Hills:
## Treatment emmean SE df lower.CL upper.CL
## Ambient 0.0 0.0817 105 -0.162 0.162
## MHW1 0.0 0.0870 105 -0.173 0.173
## MHW2 0.0 0.0870 105 -0.173 0.173
##
## Site = Cronulla:
## Treatment emmean SE df lower.CL upper.CL
## Ambient 0.0 0.0705 102 -0.140 0.140
## MHW1 0.0 0.0705 102 -0.140 0.140
## MHW2 0.0 0.0705 102 -0.140 0.140
##
## Site = Palm Beach:
## Treatment emmean SE df lower.CL upper.CL
## Ambient 0.0 0.0737 103 -0.146 0.146
## MHW1 0.0 0.0737 103 -0.146 0.146
## MHW2 0.0 0.0817 105 -0.162 0.162
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## Site = Black Head:
## contrast estimate SE df t.ratio p.value
## Ambient - MHW1 0.0 0.1131 107 0.000 1.0000
## Ambient - MHW2 -0.5 0.1127 106 -4.435 0.0001
## MHW1 - MHW2 -0.5 0.1098 107 -4.553 <.0001
##
## Site = Bonny Hills:
## contrast estimate SE df t.ratio p.value
## Ambient - MHW1 0.0 0.1191 106 0.000 1.0000
## Ambient - MHW2 0.0 0.1191 106 0.000 1.0000
## MHW1 - MHW2 0.0 0.1227 106 0.000 1.0000
##
## Site = Cronulla:
## contrast estimate SE df t.ratio p.value
## Ambient - MHW1 0.0 0.0997 104 0.000 1.0000
## Ambient - MHW2 0.0 0.0997 104 0.000 1.0000
## MHW1 - MHW2 0.0 0.0997 104 0.000 1.0000
##
## Site = Palm Beach:
## contrast estimate SE df t.ratio p.value
## Ambient - MHW1 0.0 0.1043 105 0.000 1.0000
## Ambient - MHW2 0.0 0.1099 105 0.000 1.0000
## MHW1 - MHW2 0.0 0.1099 105 0.000 1.0000
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimatesemm3 <- emmeans(my_lmer, specs = pairwise ~ Site|Treatment, type = "response", adjust = "tukey")
summary(emm3)## $emmeans
## Treatment = Ambient:
## Site emmean SE df lower.CL upper.CL
## Black Head 0.0 0.0820 103 -0.163 0.163
## Bonny Hills 0.0 0.0817 105 -0.162 0.162
## Cronulla 0.0 0.0705 102 -0.140 0.140
## Palm Beach 0.0 0.0737 103 -0.146 0.146
##
## Treatment = MHW1:
## Site emmean SE df lower.CL upper.CL
## Black Head 0.0 0.0774 104 -0.154 0.154
## Bonny Hills 0.0 0.0870 105 -0.173 0.173
## Cronulla 0.0 0.0705 102 -0.140 0.140
## Palm Beach 0.0 0.0737 103 -0.146 0.146
##
## Treatment = MHW2:
## Site emmean SE df lower.CL upper.CL
## Black Head 0.5 0.0777 102 0.346 0.654
## Bonny Hills 0.0 0.0870 105 -0.173 0.173
## Cronulla 0.0 0.0705 102 -0.140 0.140
## Palm Beach 0.0 0.0817 105 -0.162 0.162
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## Treatment = Ambient:
## contrast estimate SE df t.ratio p.value
## Black Head - Bonny Hills 0.0 0.115 105 0.000 1.0000
## Black Head - Cronulla 0.0 0.108 106 0.000 1.0000
## Black Head - Palm Beach 0.0 0.110 107 0.000 1.0000
## Bonny Hills - Cronulla 0.0 0.108 105 0.000 1.0000
## Bonny Hills - Palm Beach 0.0 0.110 105 0.000 1.0000
## Cronulla - Palm Beach 0.0 0.102 104 0.000 1.0000
##
## Treatment = MHW1:
## contrast estimate SE df t.ratio p.value
## Black Head - Bonny Hills 0.0 0.117 108 0.000 1.0000
## Black Head - Cronulla 0.0 0.105 105 0.000 1.0000
## Black Head - Palm Beach 0.0 0.107 105 0.000 1.0000
## Bonny Hills - Cronulla 0.0 0.112 106 0.000 1.0000
## Bonny Hills - Palm Beach 0.0 0.114 107 0.000 1.0000
## Cronulla - Palm Beach 0.0 0.102 104 0.000 1.0000
##
## Treatment = MHW2:
## contrast estimate SE df t.ratio p.value
## Black Head - Bonny Hills 0.5 0.116 105 4.309 0.0002
## Black Head - Cronulla 0.5 0.105 106 4.767 <.0001
## Black Head - Palm Beach 0.5 0.112 105 4.445 0.0001
## Bonny Hills - Cronulla 0.0 0.112 106 0.000 1.0000
## Bonny Hills - Palm Beach 0.0 0.119 106 0.000 1.0000
## Cronulla - Palm Beach 0.0 0.108 105 0.000 1.0000
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 4 estimates- Check model assumptions
res <- resid(my_lmer)
#produce residual vs. fitted plot
plot(fitted(my_lmer), res)
#add a horizontal line at 0
abline(0,0)
#create Q-Q plot for residuals
qqnorm(res)
#add a straight diagonal line to the plot
qqline(res) 
2.2. Day 12
- Run LMM for green colouration at day 12
my_glmer<- lmer(Green ~ Site*Treatment + (1|Location), data=Juveniles_stat12_1)## boundary (singular) fit: see help('isSingular')summary(my_glmer)## Linear mixed model fit by REML ['lmerMod']
## Formula: Green ~ Site * Treatment + (1 | Location)
## Data: Juveniles_stat12_1
##
## REML criterion at convergence: 28.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.057 0.000 0.000 0.000 6.171
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 0.00000 0.0000
## Residual 0.05909 0.2431
## Number of obs: 122, groups: Location, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -7.935e-17 8.103e-02 0.000
## SiteBonny Hills 1.490e-16 1.146e-01 0.000
## SiteCronulla 8.670e-17 1.072e-01 0.000
## SitePalm Beach 1.706e-16 1.093e-01 0.000
## TreatmentMHW1 -1.058e-17 1.093e-01 0.000
## TreatmentMHW2 5.000e-01 1.117e-01 4.477
## SiteBonny Hills:TreatmentMHW1 -5.315e-17 1.609e-01 0.000
## SiteCronulla:TreatmentMHW1 0.000e+00 1.476e-01 0.000
## SitePalm Beach:TreatmentMHW1 -1.185e-16 1.506e-01 0.000
## SiteBonny Hills:TreatmentMHW2 -5.000e-01 1.626e-01 -3.076
## SiteCronulla:TreatmentMHW2 -5.000e-01 1.494e-01 -3.347
## SitePalm Beach:TreatmentMHW2 -5.000e-01 1.562e-01 -3.200
##
## Correlation of Fixed Effects:
## (Intr) StBnnH StCrnl StPlmB TrMHW1 TrMHW2 SBH:TMHW1 SC:TMHW1
## SitBnnyHlls -0.707
## SiteCronull -0.756 0.535
## SitePalmBch -0.742 0.524 0.561
## TretmntMHW1 -0.742 0.524 0.561 0.550
## TretmntMHW2 -0.725 0.513 0.548 0.538 0.538
## StBHl:TMHW1 0.504 -0.712 -0.381 -0.373 -0.679 -0.365
## StCrn:TMHW1 0.549 -0.388 -0.726 -0.407 -0.740 -0.398 0.503
## StPBc:TMHW1 0.538 -0.380 -0.407 -0.725 -0.725 -0.390 0.493 0.537
## StBHl:TMHW2 0.498 -0.705 -0.377 -0.370 -0.370 -0.687 0.502 0.274
## StCrn:TMHW2 0.542 -0.383 -0.717 -0.402 -0.402 -0.748 0.273 0.521
## StPBc:TMHW2 0.519 -0.367 -0.392 -0.699 -0.385 -0.715 0.261 0.285
## SPB:TMHW1 SBH:TMHW2 SC:TMHW2
## SitBnnyHlls
## SiteCronull
## SitePalmBch
## TretmntMHW1
## TretmntMHW2
## StBHl:TMHW1
## StCrn:TMHW1
## StPBc:TMHW1
## StBHl:TMHW2 0.268
## StCrn:TMHW2 0.292 0.514
## StPBc:TMHW2 0.507 0.491 0.534
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')my_glmer_anova <- car::Anova(my_glmer)
my_glmer_anova## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Green
## Chisq Df Pr(>Chisq)
## Site 10.3689 3 0.015677 *
## Treatment 7.1145 2 0.028518 *
## Site:Treatment 21.0907 6 0.001767 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1overdisp.glmer(my_glmer)## Residual deviance: 6.5 on 109 degrees of freedom (ratio: 0.06)- Post-hoc test
emm3 <- emmeans(my_glmer, specs = pairwise ~ Treatment|Site, type = "response", adjust = "tukey")
summary(emm3)## $emmeans
## Site = Black Head:
## Treatment emmean SE df lower.CL upper.CL
## Ambient 0.0 0.0816 104 -0.162 0.162
## MHW1 0.0 0.0734 104 -0.146 0.146
## MHW2 0.5 0.0773 103 0.347 0.653
##
## Site = Bonny Hills:
## Treatment emmean SE df lower.CL upper.CL
## Ambient 0.0 0.0813 106 -0.161 0.161
## MHW1 0.0 0.0866 106 -0.172 0.172
## MHW2 0.0 0.0866 106 -0.172 0.172
##
## Site = Cronulla:
## Treatment emmean SE df lower.CL upper.CL
## Ambient 0.0 0.0702 103 -0.139 0.139
## MHW1 0.0 0.0702 103 -0.139 0.139
## MHW2 0.0 0.0702 103 -0.139 0.139
##
## Site = Palm Beach:
## Treatment emmean SE df lower.CL upper.CL
## Ambient 0.0 0.0734 104 -0.146 0.146
## MHW1 0.0 0.0734 104 -0.146 0.146
## MHW2 0.0 0.0813 106 -0.161 0.161
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## Site = Black Head:
## contrast estimate SE df t.ratio p.value
## Ambient - MHW1 0.0 0.1099 108 0.000 1.0000
## Ambient - MHW2 -0.5 0.1122 107 -4.456 0.0001
## MHW1 - MHW2 -0.5 0.1067 107 -4.688 <.0001
##
## Site = Bonny Hills:
## contrast estimate SE df t.ratio p.value
## Ambient - MHW1 0.0 0.1185 107 0.000 1.0000
## Ambient - MHW2 0.0 0.1185 107 0.000 1.0000
## MHW1 - MHW2 0.0 0.1221 107 0.000 1.0000
##
## Site = Cronulla:
## contrast estimate SE df t.ratio p.value
## Ambient - MHW1 0.0 0.0992 105 0.000 1.0000
## Ambient - MHW2 0.0 0.0992 105 0.000 1.0000
## MHW1 - MHW2 0.0 0.0992 105 0.000 1.0000
##
## Site = Palm Beach:
## contrast estimate SE df t.ratio p.value
## Ambient - MHW1 0.0 0.1038 106 0.000 1.0000
## Ambient - MHW2 0.0 0.1094 106 0.000 1.0000
## MHW1 - MHW2 0.0 0.1094 106 0.000 1.0000
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimatesemm3 <- emmeans(my_glmer, specs = pairwise ~ Site|Treatment, type = "response", adjust = "tukey")
summary(emm3)## $emmeans
## Treatment = Ambient:
## Site emmean SE df lower.CL upper.CL
## Black Head 0.0 0.0816 104 -0.162 0.162
## Bonny Hills 0.0 0.0813 106 -0.161 0.161
## Cronulla 0.0 0.0702 103 -0.139 0.139
## Palm Beach 0.0 0.0734 104 -0.146 0.146
##
## Treatment = MHW1:
## Site emmean SE df lower.CL upper.CL
## Black Head 0.0 0.0734 104 -0.146 0.146
## Bonny Hills 0.0 0.0866 106 -0.172 0.172
## Cronulla 0.0 0.0702 103 -0.139 0.139
## Palm Beach 0.0 0.0734 104 -0.146 0.146
##
## Treatment = MHW2:
## Site emmean SE df lower.CL upper.CL
## Black Head 0.5 0.0773 103 0.347 0.653
## Bonny Hills 0.0 0.0866 106 -0.172 0.172
## Cronulla 0.0 0.0702 103 -0.139 0.139
## Palm Beach 0.0 0.0813 106 -0.161 0.161
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## Treatment = Ambient:
## contrast estimate SE df t.ratio p.value
## Black Head - Bonny Hills 0.0 0.115 106 0.000 1.0000
## Black Head - Cronulla 0.0 0.108 107 0.000 1.0000
## Black Head - Palm Beach 0.0 0.110 108 0.000 1.0000
## Bonny Hills - Cronulla 0.0 0.107 106 0.000 1.0000
## Bonny Hills - Palm Beach 0.0 0.109 106 0.000 1.0000
## Cronulla - Palm Beach 0.0 0.102 106 0.000 1.0000
##
## Treatment = MHW1:
## contrast estimate SE df t.ratio p.value
## Black Head - Bonny Hills 0.0 0.114 108 0.000 1.0000
## Black Head - Cronulla 0.0 0.102 106 0.000 1.0000
## Black Head - Palm Beach 0.0 0.104 106 0.000 1.0000
## Bonny Hills - Cronulla 0.0 0.111 107 0.000 1.0000
## Bonny Hills - Palm Beach 0.0 0.114 108 0.000 1.0000
## Cronulla - Palm Beach 0.0 0.102 106 0.000 1.0000
##
## Treatment = MHW2:
## contrast estimate SE df t.ratio p.value
## Black Head - Bonny Hills 0.5 0.116 106 4.329 0.0002
## Black Head - Cronulla 0.5 0.104 107 4.789 <.0001
## Black Head - Palm Beach 0.5 0.112 106 4.465 0.0001
## Bonny Hills - Cronulla 0.0 0.111 107 0.000 1.0000
## Bonny Hills - Palm Beach 0.0 0.119 107 0.000 1.0000
## Cronulla - Palm Beach 0.0 0.107 106 0.000 1.0000
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 4 estimates- Check model assumptions
res <- resid(my_glmer)
#produce residual vs. fitted plot
plot(fitted(my_glmer), res)
#add a horizontal line at 0
abline(0,0)
#create Q-Q plot for residuals
qqnorm(res)
#add a straight diagonal line to the plot
qqline(res) 
2.3. Day 18
- Run LMM for green colouration at day 18
my_glmer1<- lmer(Green ~ Site*Treatment + (1|Location), data=Juveniles_stat18_1)## boundary (singular) fit: see help('isSingular')summary(my_glmer1)## Linear mixed model fit by REML ['lmerMod']
## Formula: Green ~ Site * Treatment + (1 | Location)
## Data: Juveniles_stat18_1
##
## REML criterion at convergence: 193.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.1968 -0.4660 -0.1898 -0.1424 2.9512
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 0.0000 0.0000
## Residual 0.3425 0.5853
## Number of obs: 107, groups: Location, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.22222 0.19509 1.139
## SiteBonny Hills -0.11111 0.27590 -0.403
## SiteCronulla -0.12222 0.26891 -0.455
## SitePalm Beach 0.05051 0.26306 0.192
## TreatmentMHW1 0.07778 0.26891 0.289
## TreatmentMHW2 0.27778 0.30846 0.901
## SiteBonny Hills:TreatmentMHW1 -0.04603 0.39913 -0.115
## SiteCronulla:TreatmentMHW1 -0.09444 0.36758 -0.257
## SitePalm Beach:TreatmentMHW1 -0.25960 0.36687 -0.708
## SiteBonny Hills:TreatmentMHW2 0.89683 0.42678 2.101
## SiteCronulla:TreatmentMHW2 0.24722 0.41500 0.596
## SitePalm Beach:TreatmentMHW2 -0.26479 0.41860 -0.633
##
## Correlation of Fixed Effects:
## (Intr) StBnnH StCrnl StPlmB TrMHW1 TrMHW2 SBH:TMHW1 SC:TMHW1
## SitBnnyHlls -0.707
## SiteCronull -0.725 0.513
## SitePalmBch -0.742 0.524 0.538
## TretmntMHW1 -0.725 0.513 0.526 0.538
## TretmntMHW2 -0.632 0.447 0.459 0.469 0.459
## StBHl:TMHW1 0.489 -0.691 -0.355 -0.362 -0.674 -0.309
## StCrn:TMHW1 0.531 -0.375 -0.732 -0.394 -0.732 -0.336 0.493
## StPBc:TMHW1 0.532 -0.376 -0.386 -0.717 -0.733 -0.336 0.494 0.536
## StBHl:TMHW2 0.457 -0.646 -0.332 -0.339 -0.332 -0.723 0.447 0.243
## StCrn:TMHW2 0.470 -0.332 -0.648 -0.349 -0.341 -0.743 0.230 0.474
## StPBc:TMHW2 0.466 -0.330 -0.338 -0.628 -0.338 -0.737 0.228 0.247
## SPB:TMHW1 SBH:TMHW2 SC:TMHW2
## SitBnnyHlls
## SiteCronull
## SitePalmBch
## TretmntMHW1
## TretmntMHW2
## StBHl:TMHW1
## StCrn:TMHW1
## StPBc:TMHW1
## StBHl:TMHW2 0.243
## StCrn:TMHW2 0.250 0.537
## StPBc:TMHW2 0.451 0.533 0.548
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')overdisp.glmer(my_glmer1)## Residual deviance: 32.541 on 94 degrees of freedom (ratio: 0.346)my_glmer_anova <- car::Anova(my_glmer1)
my_glmer_anova## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Green
## Chisq Df Pr(>Chisq)
## Site 2.6847 3 0.442830
## Treatment 15.3947 2 0.000454 ***
## Site:Treatment 10.1880 6 0.116955
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- Post-hoc test
emm3 <- emmeans(my_glmer1, specs = pairwise ~ Treatment, type = "response", adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionssummary(emm3)## $emmeans
## Treatment emmean SE df lower.CL upper.CL
## Ambient 0.177 0.0949 28.9 -0.0177 0.371
## MHW1 0.154 0.0950 30.6 -0.0396 0.348
## MHW2 0.674 0.1138 36.8 0.4435 0.905
##
## Results are averaged over the levels of: Site
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## Ambient - MHW1 0.0222 0.134 92.2 0.166 0.9849
## Ambient - MHW2 -0.4976 0.147 93.6 -3.391 0.0029
## MHW1 - MHW2 -0.5198 0.147 94.1 -3.531 0.0018
##
## Results are averaged over the levels of: Site
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimates- Check model assumptions
res <- resid(my_glmer1)
#produce residual vs. fitted plot
plot(fitted(my_glmer1), res)
#add a horizontal line at 0
abline(0,0)
#create Q-Q plot for residuals
qqnorm(res)
#add a straight diagonal line to the plot
qqline(res) 
2.4. Day 22
- Run LMM for green colouration at day 22
my_glmer1<- lmer(Green~ Site*Treatment + (1|Location), data=Juveniles_stat22_1)## boundary (singular) fit: see help('isSingular')summary(my_glmer1)## Linear mixed model fit by REML ['lmerMod']
## Formula: Green ~ Site * Treatment + (1 | Location)
## Data: Juveniles_stat22_1
##
## REML criterion at convergence: 188.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.0666 -0.5333 -0.2400 0.1333 3.1198
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 0.0000 0.0000
## Residual 0.6945 0.8334
## Number of obs: 77, groups: Location, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.44444 0.27780 1.600
## SiteBonny Hills 0.22222 0.39287 0.566
## SiteCronulla -0.04444 0.38292 -0.116
## SitePalm Beach -0.17172 0.37458 -0.458
## TreatmentMHW1 0.44444 0.39287 1.131
## SiteBonny Hills:TreatmentMHW1 -0.53968 0.57510 -0.938
## SiteCronulla:TreatmentMHW1 -0.76111 0.53073 -1.434
## SitePalm Beach:TreatmentMHW1 -0.51717 0.53567 -0.965
##
## Correlation of Fixed Effects:
## (Intr) StBnnH StCrnl StPlmB TrMHW1 SBH:TM SC:TMH
## SitBnnyHlls -0.707
## SiteCronull -0.725 0.513
## SitePalmBch -0.742 0.524 0.538
## TretmntMHW1 -0.707 0.500 0.513 0.524
## StBHl:TMHW1 0.483 -0.683 -0.350 -0.358 -0.683
## StCrn:TMHW1 0.523 -0.370 -0.721 -0.388 -0.740 0.506
## StPBc:TMHW1 0.519 -0.367 -0.376 -0.699 -0.733 0.501 0.543
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')overdisp.glmer(my_glmer1)## Residual deviance: 47.924 on 68 degrees of freedom (ratio: 0.705)my_glmer_anova <- car::Anova(my_glmer1)
my_glmer_anova## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Green
## Chisq Df Pr(>Chisq)
## Site 4.6904 3 0.1959
## Treatment 0.0176 1 0.8945
## Site:Treatment 2.1411 3 0.5437- Check model assumptions
res <- resid(my_glmer1)
#produce residual vs. fitted plot
plot(fitted(my_glmer1), res)
#add a horizontal line at 0
abline(0,0)
#create Q-Q plot for residuals
qqnorm(res)
#add a straight diagonal line to the plot
qqline(res) 
3f) Necrosis
1. Visualise differences
Necrosis_plot <- Juv_means %>%
filter(Day_equivalent == 8|Day_equivalent == 12|Day_equivalent == 18| Day_equivalent == 22) %>%
ggplot(aes(x= Day_equivalent, y= Mean_necrosis, colour = Site)) +
geom_errorbar(aes(ymin=Mean_necrosis-se_necrosis, ymax=Mean_necrosis+se_necrosis),colour = "grey", width=.2) +
geom_point() +
geom_line()+
Theme1 +
labs( y="Necrosis (0-4)", x="Time (days)") +
scale_colour_manual(values = c("darkred", "orange", "darkolivegreen3", "royalblue3")) +
theme(strip.text.x = element_text(size = 12))+
scale_y_continuous(limits = c(0, 4.5), breaks = seq(0, 4, by = 1))+
scale_x_continuous(limits = c(7, 23), breaks = c(8, 12, 18,22))+
facet_wrap(~ factor(Treatment, c("Ambient","MHW1","MHW2")), ncol = 3)
Necrosis_plot
ggsave("Necrosis_juv.jpeg", width = 220, height = 60, units = "mm")2. Statistical analyses
2.1. Day 8
- Run LMM for necrosis at day 8
my_lmer<- lmer(Black ~ Site*Treatment + (1|Location), data=Juveniles_stat8_1)
summary(my_lmer)## Linear mixed model fit by REML ['lmerMod']
## Formula: Black ~ Site * Treatment + (1 | Location)
## Data: Juveniles_stat8_1
##
## REML criterion at convergence: 170.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.2562 -0.3858 -0.1838 0.0310 3.6470
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 0.00869 0.09322
## Residual 0.21083 0.45916
## Number of obs: 121, groups: Location, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.102387 0.158384 0.646
## SiteBonny Hills -0.095433 0.216736 -0.440
## SiteCronulla -0.102387 0.202993 -0.504
## SitePalm Beach 0.084893 0.207102 0.410
## TreatmentMHW1 -0.008753 0.211872 -0.041
## TreatmentMHW2 0.414317 0.211615 1.958
## SiteBonny Hills:TreatmentMHW1 0.152393 0.307582 0.495
## SiteCronulla:TreatmentMHW1 0.092086 0.282892 0.326
## SitePalm Beach:TreatmentMHW1 0.086835 0.288972 0.300
## SiteBonny Hills:TreatmentMHW2 -0.395530 0.307182 -1.288
## SiteCronulla:TreatmentMHW2 0.085683 0.282700 0.303
## SitePalm Beach:TreatmentMHW2 -0.483532 0.295861 -1.634
##
## Correlation of Fixed Effects:
## (Intr) StBnnH StCrnl StPlmB TrMHW1 TrMHW2 SBH:TMHW1 SC:TMHW1
## SitBnnyHlls -0.686
## SiteCronull -0.735 0.535
## SitePalmBch -0.721 0.526 0.563
## TretmntMHW1 -0.706 0.513 0.551 0.540
## TretmntMHW2 -0.704 0.513 0.549 0.538 0.526
## StBHl:TMHW1 0.483 -0.703 -0.377 -0.369 -0.687 -0.360
## StCrn:TMHW1 0.529 -0.384 -0.718 -0.405 -0.749 -0.394 0.514
## StPBc:TMHW1 0.520 -0.378 -0.406 -0.719 -0.735 -0.388 0.503 0.550
## StBHl:TMHW2 0.483 -0.705 -0.377 -0.370 -0.361 -0.686 0.496 0.270
## StCrn:TMHW2 0.527 -0.384 -0.717 -0.403 -0.394 -0.749 0.270 0.515
## StPBc:TMHW2 0.506 -0.368 -0.395 -0.700 -0.379 -0.716 0.259 0.284
## SPB:TMHW1 SBH:TMHW2 SC:TMHW2
## SitBnnyHlls
## SiteCronull
## SitePalmBch
## TretmntMHW1
## TretmntMHW2
## StBHl:TMHW1
## StCrn:TMHW1
## StPBc:TMHW1
## StBHl:TMHW2 0.265
## StCrn:TMHW2 0.290 0.513
## StPBc:TMHW2 0.504 0.491 0.536my_glmer_anova <- car::Anova(my_lmer)
my_glmer_anova## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Black
## Chisq Df Pr(>Chisq)
## Site 2.1995 3 0.53205
## Treatment 5.7962 2 0.05513 .
## Site:Treatment 8.7888 6 0.18581
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- Check model assumptions
res <- resid(my_lmer)
#produce residual vs. fitted plot
plot(fitted(my_lmer), res)
#add a horizontal line at 0
abline(0,0)
#create Q-Q plot for residuals
qqnorm(res)
#add a straight diagonal line to the plot
qqline(res) 
2.2. Day 12
- Run LMM for necrosis at day 12
my_glmer<- lmer(Black ~ Site*Treatment + (1|Location), data=Juveniles_stat12_1)
summary(my_glmer)## Linear mixed model fit by REML ['lmerMod']
## Formula: Black ~ Site * Treatment + (1 | Location)
## Data: Juveniles_stat12_1
##
## REML criterion at convergence: 186.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.2525 -0.4755 -0.1819 0.0069 3.4461
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 0.002206 0.04697
## Residual 0.246803 0.49679
## Number of obs: 122, groups: Location, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.108318 0.166949 0.649
## SiteBonny Hills -0.106455 0.234291 -0.454
## SiteCronulla -0.108318 0.219251 -0.494
## SitePalm Beach 0.075089 0.223548 0.336
## TreatmentMHW1 -0.018640 0.223546 -0.083
## TreatmentMHW2 0.496421 0.228487 2.173
## SiteBonny Hills:TreatmentMHW1 0.398622 0.328957 1.212
## SiteCronulla:TreatmentMHW1 0.101973 0.301839 0.338
## SitePalm Beach:TreatmentMHW1 0.105740 0.308146 0.343
## SiteBonny Hills:TreatmentMHW2 -0.491866 0.332273 -1.480
## SiteCronulla:TreatmentMHW2 0.003579 0.305516 0.012
## SitePalm Beach:TreatmentMHW2 -0.455743 0.319577 -1.426
##
## Correlation of Fixed Effects:
## (Intr) StBnnH StCrnl StPlmB TrMHW1 TrMHW2 SBH:TMHW1 SC:TMHW1
## SitBnnyHlls -0.702
## SiteCronull -0.751 0.535
## SitePalmBch -0.737 0.525 0.561
## TretmntMHW1 -0.737 0.524 0.561 0.551
## TretmntMHW2 -0.720 0.513 0.549 0.538 0.538
## StBHl:TMHW1 0.500 -0.712 -0.381 -0.373 -0.679 -0.365
## StCrn:TMHW1 0.546 -0.388 -0.727 -0.408 -0.741 -0.399 0.503
## StPBc:TMHW1 0.536 -0.381 -0.408 -0.726 -0.726 -0.391 0.492 0.538
## StBHl:TMHW2 0.495 -0.705 -0.377 -0.370 -0.369 -0.687 0.502 0.274
## StCrn:TMHW2 0.539 -0.384 -0.717 -0.402 -0.402 -0.748 0.273 0.521
## StPBc:TMHW2 0.516 -0.367 -0.393 -0.699 -0.385 -0.715 0.261 0.285
## SPB:TMHW1 SBH:TMHW2 SC:TMHW2
## SitBnnyHlls
## SiteCronull
## SitePalmBch
## TretmntMHW1
## TretmntMHW2
## StBHl:TMHW1
## StCrn:TMHW1
## StPBc:TMHW1
## StBHl:TMHW2 0.268
## StCrn:TMHW2 0.292 0.514
## StPBc:TMHW2 0.508 0.491 0.535my_glmer_anova <- car::Anova(my_glmer)
my_glmer_anova## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Black
## Chisq Df Pr(>Chisq)
## Site 1.1086 3 0.77499
## Treatment 6.7140 2 0.03484 *
## Site:Treatment 10.5056 6 0.10491
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1overdisp.glmer(my_glmer)## Residual deviance: 26.962 on 109 degrees of freedom (ratio: 0.247)- Post-hoc test
emm3 <- emmeans(my_glmer, specs = pairwise ~ Treatment, type = "response", adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionssummary(emm3)## $emmeans
## Treatment emmean SE df lower.CL upper.CL
## Ambient 0.0734 0.0810 30.2 -0.0919 0.239
## MHW1 0.2063 0.0800 30.4 0.0429 0.370
## MHW2 0.3338 0.0835 31.4 0.1636 0.504
##
## Results are averaged over the levels of: Site
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## Ambient - MHW1 -0.133 0.110 106 -1.204 0.4535
## Ambient - MHW2 -0.260 0.113 107 -2.312 0.0585
## MHW1 - MHW2 -0.127 0.112 107 -1.135 0.4945
##
## Results are averaged over the levels of: Site
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimatesCheck model assumptions
res <- resid(my_glmer)
#produce residual vs. fitted plot
plot(fitted(my_glmer), res)
#add a horizontal line at 0
abline(0,0)
#create Q-Q plot for residuals
qqnorm(res)
#add a straight diagonal line to the plot
qqline(res) 
2.3. Day 18
- Run LMM for necrosis at day 18
my_glmer1<- lmer(Black ~ Site*Treatment + (1|Location), data=Juveniles_stat18_1)## boundary (singular) fit: see help('isSingular')summary(my_glmer1)## Linear mixed model fit by REML ['lmerMod']
## Formula: Black ~ Site * Treatment + (1 | Location)
## Data: Juveniles_stat18_1
##
## REML criterion at convergence: 213.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.3988 -0.6839 -0.2198 0.6924 3.6928
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 0.0000 0.0000
## Residual 0.4224 0.6499
## Number of obs: 107, groups: Location, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.44444 0.21664 2.052
## SiteBonny Hills -0.11111 0.30637 -0.363
## SiteCronulla -0.14444 0.29861 -0.484
## SitePalm Beach 0.46465 0.29211 1.591
## TreatmentMHW1 0.15556 0.29861 0.521
## TreatmentMHW2 0.05556 0.34253 0.162
## SiteBonny Hills:TreatmentMHW1 0.22540 0.44322 0.509
## SiteCronulla:TreatmentMHW1 -0.20556 0.40818 -0.504
## SitePalm Beach:TreatmentMHW1 -0.61010 0.40739 -1.498
## SiteBonny Hills:TreatmentMHW2 -0.24603 0.47392 -0.519
## SiteCronulla:TreatmentMHW2 -0.23056 0.46083 -0.500
## SitePalm Beach:TreatmentMHW2 -0.10750 0.46483 -0.231
##
## Correlation of Fixed Effects:
## (Intr) StBnnH StCrnl StPlmB TrMHW1 TrMHW2 SBH:TMHW1 SC:TMHW1
## SitBnnyHlls -0.707
## SiteCronull -0.725 0.513
## SitePalmBch -0.742 0.524 0.538
## TretmntMHW1 -0.725 0.513 0.526 0.538
## TretmntMHW2 -0.632 0.447 0.459 0.469 0.459
## StBHl:TMHW1 0.489 -0.691 -0.355 -0.362 -0.674 -0.309
## StCrn:TMHW1 0.531 -0.375 -0.732 -0.394 -0.732 -0.336 0.493
## StPBc:TMHW1 0.532 -0.376 -0.386 -0.717 -0.733 -0.336 0.494 0.536
## StBHl:TMHW2 0.457 -0.646 -0.332 -0.339 -0.332 -0.723 0.447 0.243
## StCrn:TMHW2 0.470 -0.332 -0.648 -0.349 -0.341 -0.743 0.230 0.474
## StPBc:TMHW2 0.466 -0.330 -0.338 -0.628 -0.338 -0.737 0.228 0.247
## SPB:TMHW1 SBH:TMHW2 SC:TMHW2
## SitBnnyHlls
## SiteCronull
## SitePalmBch
## TretmntMHW1
## TretmntMHW2
## StBHl:TMHW1
## StCrn:TMHW1
## StPBc:TMHW1
## StBHl:TMHW2 0.243
## StCrn:TMHW2 0.250 0.537
## StPBc:TMHW2 0.451 0.533 0.548
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')overdisp.glmer(my_glmer1)## Residual deviance: 40.126 on 94 degrees of freedom (ratio: 0.427)my_glmer_anova <- car::Anova(my_glmer1)
my_glmer_anova## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Black
## Chisq Df Pr(>Chisq)
## Site 8.7464 3 0.03286 *
## Treatment 0.4148 2 0.81271
## Site:Treatment 6.0997 6 0.41212
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- Post-hoc test
emm3 <- emmeans(my_glmer1, specs = pairwise ~ Site, type = "response", adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionssummary(emm3)## $emmeans
## Site emmean SE df lower.CL upper.CL
## Black Head 0.515 0.134 54.6 0.246 0.784
## Bonny Hills 0.397 0.139 50.1 0.118 0.676
## Cronulla 0.225 0.121 43.4 -0.019 0.469
## Palm Beach 0.740 0.124 47.5 0.491 0.990
##
## Results are averaged over the levels of: Treatment
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## Black Head - Bonny Hills 0.118 0.191 92.1 0.616 0.9266
## Black Head - Cronulla 0.290 0.180 91.2 1.611 0.3776
## Black Head - Palm Beach -0.225 0.182 91.0 -1.239 0.6038
## Bonny Hills - Cronulla 0.172 0.183 93.1 0.938 0.7846
## Bonny Hills - Palm Beach -0.343 0.185 92.1 -1.859 0.2529
## Cronulla - Palm Beach -0.515 0.173 91.1 -2.984 0.0188
##
## Results are averaged over the levels of: Treatment
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 4 estimates- Check model assumptions
res <- resid(my_glmer1)
#produce residual vs. fitted plot
plot(fitted(my_glmer1), res)
#add a horizontal line at 0
abline(0,0)
#create Q-Q plot for residuals
qqnorm(res)
#add a straight diagonal line to the plot
qqline(res) 
2.4. Day 22
- Run LMM for necrosis at day 22
my_glmer1<- lmer(Black~ Site*Treatment + (1|Location), data=Juveniles_stat22_1)## boundary (singular) fit: see help('isSingular')summary(my_glmer1)## Linear mixed model fit by REML ['lmerMod']
## Formula: Black ~ Site * Treatment + (1 | Location)
## Data: Juveniles_stat22_1
##
## REML criterion at convergence: 207.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.4962 -0.6284 -0.1904 0.4655 2.6932
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 8.337e-18 2.887e-09
## Residual 9.116e-01 9.548e-01
## Number of obs: 77, groups: Location, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 5.556e-01 3.183e-01 1.746
## SiteBonny Hills -1.091e-16 4.501e-01 0.000
## SiteCronulla 4.444e-02 4.387e-01 0.101
## SitePalm Beach 6.263e-01 4.291e-01 1.459
## TreatmentMHW1 2.222e-01 4.501e-01 0.494
## SiteBonny Hills:TreatmentMHW1 6.508e-01 6.589e-01 0.988
## SiteCronulla:TreatmentMHW1 -7.222e-02 6.080e-01 -0.119
## SitePalm Beach:TreatmentMHW1 -3.040e-01 6.137e-01 -0.495
##
## Correlation of Fixed Effects:
## (Intr) StBnnH StCrnl StPlmB TrMHW1 SBH:TM SC:TMH
## SitBnnyHlls -0.707
## SiteCronull -0.725 0.513
## SitePalmBch -0.742 0.524 0.538
## TretmntMHW1 -0.707 0.500 0.513 0.524
## StBHl:TMHW1 0.483 -0.683 -0.350 -0.358 -0.683
## StCrn:TMHW1 0.523 -0.370 -0.721 -0.388 -0.740 0.506
## StPBc:TMHW1 0.519 -0.367 -0.376 -0.699 -0.733 0.501 0.543
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')overdisp.glmer(my_glmer1)## Residual deviance: 62.901 on 68 degrees of freedom (ratio: 0.925)my_glmer_anova <- car::Anova(my_glmer1)
my_glmer_anova## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Black
## Chisq Df Pr(>Chisq)
## Site 3.6704 3 0.2993
## Treatment 1.3361 1 0.2477
## Site:Treatment 2.3727 3 0.4987- Check model assumptions
res <- resid(my_glmer1)
#produce residual vs. fitted plot
plot(fitted(my_glmer1), res)
#add a horizontal line at 0
abline(0,0)
#create Q-Q plot for residuals
qqnorm(res)
#add a straight diagonal line to the plot
qqline(res) 
3g) Photosynthetic efficiency
1. Visualise differences
FvFm_plot <- Photo_means %>%
ggplot(aes(x= Day_equivalent, y= Mean_FvFm, colour = Site)) +
geom_errorbar(aes(ymin=Mean_FvFm-se_FvFm, ymax=Mean_FvFm+se_FvFm),colour = "grey", width=.2) +
geom_point() +
geom_line()+
Theme1 +
labs( y="Fv/Fm", x="Time (days)") +
scale_colour_manual(values = c("darkred", "orange", "darkolivegreen3", "royalblue3")) +
theme(strip.text.x = element_text(size = 12))+
scale_x_continuous(limits = c(7, 23), breaks = c(8, 12, 18,22))+
facet_wrap(~ factor(Treatment, c("Ambient","MHW1","MHW2")), ncol = 3)
Photo_means$Day_equivalent <- as.numeric(Photo_means$Day_equivalent)
FvFm_plot <- Photo_means %>%
ggplot(aes(x= Day_equivalent, y=Mean_FvFm, fill = Treatment)) +
geom_bar(stat="identity", color="black",
position=position_dodge()) +
geom_errorbar(aes(ymin=Mean_FvFm-se_FvFm, ymax=Mean_FvFm+se_FvFm), width=.2,
position=position_dodge(0.9))+
Theme1 +
labs( y="Fv/Fm", x="Day") +
facet_wrap(~ factor(Site, c("Bonny Hills","Black Head","Palm Beach","Cronulla")), ncol = 4) +
scale_fill_manual(breaks = c("Ambient", "MHW1", "MHW2"),
values=c("yellowgreen", "tan2", "plum4"))
FvFm_plot
ggsave("photosynthesis_juv.jpeg", width = 220, height = 60, units = "mm")2. Statistical analyses
2.1. Day 8
- Run LMM for photosynthetic efficiency at day 8
my_lmer<- lmer(Fv_Fm ~ Site*Treatment + (1|Location), data=Juveniles_stat8_1)
summary(my_lmer)## Linear mixed model fit by REML ['lmerMod']
## Formula: Fv_Fm ~ Site * Treatment + (1 | Location)
## Data: Juveniles_stat8_1
##
## REML criterion at convergence: -97
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.9255 -0.6606 -0.1445 0.5488 2.9051
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 0.0008462 0.02909
## Residual 0.0177656 0.13329
## Number of obs: 119, groups: Location, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.3223903 0.0488391 6.601
## SiteBonny Hills 0.0724598 0.0649163 1.116
## SiteCronulla 0.0105264 0.0610303 0.172
## SitePalm Beach 0.0002837 0.0622072 0.005
## TreatmentMHW1 0.0032988 0.0634863 0.052
## TreatmentMHW2 0.0414982 0.0634764 0.654
## SiteBonny Hills:TreatmentMHW1 -0.0626829 0.0906227 -0.692
## SiteCronulla:TreatmentMHW1 0.0384512 0.0836149 0.460
## SitePalm Beach:TreatmentMHW1 -0.0071139 0.0865142 -0.082
## SiteBonny Hills:TreatmentMHW2 -0.1602109 0.0905539 -1.769
## SiteCronulla:TreatmentMHW2 -0.0642482 0.0836073 -0.768
## SitePalm Beach:TreatmentMHW2 0.0267890 0.0873586 0.307
##
## Correlation of Fixed Effects:
## (Intr) StBnnH StCrnl StPlmB TrMHW1 TrMHW2 SBH:TMHW1 SC:TMHW1
## SitBnnyHlls -0.706
## SiteCronull -0.753 0.565
## SitePalmBch -0.740 0.556 0.592
## TretmntMHW1 -0.724 0.542 0.580 0.569
## TretmntMHW2 -0.723 0.543 0.578 0.568 0.556
## StBHl:TMHW1 0.503 -0.713 -0.403 -0.395 -0.698 -0.386
## StCrn:TMHW1 0.550 -0.412 -0.730 -0.432 -0.759 -0.422 0.530
## StPBc:TMHW1 0.536 -0.401 -0.429 -0.722 -0.738 -0.412 0.512 0.560
## StBHl:TMHW2 0.504 -0.715 -0.403 -0.396 -0.387 -0.697 0.510 0.294
## StCrn:TMHW2 0.549 -0.413 -0.729 -0.431 -0.422 -0.759 0.293 0.532
## StPBc:TMHW2 0.527 -0.396 -0.422 -0.712 -0.406 -0.727 0.282 0.308
## SPB:TMHW1 SBH:TMHW2 SC:TMHW2
## SitBnnyHlls
## SiteCronull
## SitePalmBch
## TretmntMHW1
## TretmntMHW2
## StBHl:TMHW1
## StCrn:TMHW1
## StPBc:TMHW1
## StBHl:TMHW2 0.285
## StCrn:TMHW2 0.312 0.529
## StPBc:TMHW2 0.515 0.507 0.552my_glmer_anova <- car::Anova(my_lmer)
my_glmer_anova## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Fv_Fm
## Chisq Df Pr(>Chisq)
## Site 0.0136 3 0.9996
## Treatment 0.0481 2 0.9762
## Site:Treatment 7.0699 6 0.3144- Check model assumptions
res <- resid(my_lmer)
#produce residual vs. fitted plot
plot(fitted(my_lmer), res)
#add a horizontal line at 0
abline(0,0)
#create Q-Q plot for residuals
qqnorm(res)
#add a straight diagonal line to the plot
qqline(res) 
2.2. Day 12
- Run LMM for photosynthetic efficiency at day 12
my_glmer<- lmer(Fv_Fm ~ Site*Treatment + (1|Location), data=Juveniles_stat12_1)
summary(my_glmer)## Linear mixed model fit by REML ['lmerMod']
## Formula: Fv_Fm ~ Site * Treatment + (1 | Location)
## Data: Juveniles_stat12_1
##
## REML criterion at convergence: -105.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.88331 -0.67063 -0.08341 0.67070 3.09137
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 0.009179 0.09581
## Residual 0.015165 0.12315
## Number of obs: 120, groups: Location, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.395555 0.056995 6.940
## SiteBonny Hills 0.029451 0.058227 0.506
## SiteCronulla 0.026195 0.054610 0.480
## SitePalm Beach 0.011625 0.055787 0.208
## TreatmentMHW1 -0.007157 0.055760 -0.128
## TreatmentMHW2 0.046812 0.058582 0.799
## SiteBonny Hills:TreatmentMHW1 -0.024015 0.081695 -0.294
## SiteCronulla:TreatmentMHW1 -0.030677 0.075078 -0.409
## SitePalm Beach:TreatmentMHW1 -0.030540 0.078139 -0.391
## SiteBonny Hills:TreatmentMHW2 -0.094983 0.083525 -1.137
## SiteCronulla:TreatmentMHW2 -0.126645 0.077197 -1.641
## SitePalm Beach:TreatmentMHW2 -0.059208 0.080761 -0.733
##
## Correlation of Fixed Effects:
## (Intr) StBnnH StCrnl StPlmB TrMHW1 TrMHW2 SBH:TMHW1 SC:TMHW1
## SitBnnyHlls -0.514
## SiteCronull -0.552 0.536
## SitePalmBch -0.542 0.528 0.566
## TretmntMHW1 -0.542 0.525 0.566 0.555
## TretmntMHW2 -0.514 0.503 0.536 0.527 0.525
## StBHl:TMHW1 0.364 -0.708 -0.380 -0.371 -0.677 -0.351
## StCrn:TMHW1 0.402 -0.390 -0.728 -0.412 -0.743 -0.390 0.503
## StPBc:TMHW1 0.394 -0.381 -0.411 -0.721 -0.721 -0.384 0.483 0.535
## StBHl:TMHW2 0.357 -0.697 -0.373 -0.367 -0.364 -0.694 0.492 0.270
## StCrn:TMHW2 0.390 -0.382 -0.707 -0.400 -0.398 -0.759 0.266 0.514
## StPBc:TMHW2 0.376 -0.367 -0.393 -0.693 -0.385 -0.727 0.256 0.286
## SPB:TMHW1 SBH:TMHW2 SC:TMHW2
## SitBnnyHlls
## SiteCronull
## SitePalmBch
## TretmntMHW1
## TretmntMHW2
## StBHl:TMHW1
## StCrn:TMHW1
## StPBc:TMHW1
## StBHl:TMHW2 0.263
## StCrn:TMHW2 0.292 0.527
## StPBc:TMHW2 0.503 0.504 0.552my_glmer_anova <- car::Anova(my_glmer)
my_glmer_anova## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Fv_Fm
## Chisq Df Pr(>Chisq)
## Site 0.7408 3 0.8636
## Treatment 1.4825 2 0.4765
## Site:Treatment 3.2675 6 0.7746overdisp.glmer(my_glmer)## Residual deviance: 1.568 on 107 degrees of freedom (ratio: 0.015)- Check model assumptions
res <- resid(my_glmer)
#produce residual vs. fitted plot
plot(fitted(my_glmer), res)
#add a horizontal line at 0
abline(0,0)
#create Q-Q plot for residuals
qqnorm(res)
#add a straight diagonal line to the plot
qqline(res) 
2.3. Day 18
- Run LMM for photosynthetic efficiency at day 18
my_glmer1<- lmer(Fv_Fm ~ Site*Treatment + (1|Location), data=Juveniles_stat18_1)## boundary (singular) fit: see help('isSingular')summary(my_glmer1)## Linear mixed model fit by REML ['lmerMod']
## Formula: Fv_Fm ~ Site * Treatment + (1 | Location)
## Data: Juveniles_stat18_1
##
## REML criterion at convergence: -134.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.67047 -0.66720 0.01309 0.51564 2.71723
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 0.00000 0.0000
## Residual 0.01218 0.1104
## Number of obs: 115, groups: Location, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.257750 0.039020 6.606
## SiteBonny Hills 0.066806 0.053628 1.246
## SiteCronulla 0.103361 0.053628 1.927
## SitePalm Beach 0.036977 0.051282 0.721
## TreatmentMHW1 0.020705 0.051282 0.404
## TreatmentMHW2 -0.170639 0.053628 -3.182
## SiteBonny Hills:TreatmentMHW1 -0.084760 0.074201 -1.142
## SiteCronulla:TreatmentMHW1 -0.003399 0.070699 -0.048
## SitePalm Beach:TreatmentMHW1 -0.023795 0.069603 -0.342
## SiteBonny Hills:TreatmentMHW2 -0.059631 0.077262 -0.772
## SiteCronulla:TreatmentMHW2 -0.136199 0.073053 -1.864
## SitePalm Beach:TreatmentMHW2 0.045578 0.073053 0.624
##
## Correlation of Fixed Effects:
## (Intr) StBnnH StCrnl StPlmB TrMHW1 TrMHW2 SBH:TMHW1 SC:TMHW1
## SitBnnyHlls -0.728
## SiteCronull -0.728 0.529
## SitePalmBch -0.761 0.554 0.554
## TretmntMHW1 -0.761 0.554 0.554 0.579
## TretmntMHW2 -0.728 0.529 0.529 0.554 0.554
## StBHl:TMHW1 0.526 -0.723 -0.383 -0.400 -0.691 -0.383
## StCrn:TMHW1 0.552 -0.402 -0.759 -0.420 -0.725 -0.402 0.501
## StPBc:TMHW1 0.561 -0.408 -0.408 -0.737 -0.737 -0.408 0.509 0.534
## StBHl:TMHW2 0.505 -0.694 -0.367 -0.384 -0.384 -0.694 0.502 0.279
## StCrn:TMHW2 0.534 -0.389 -0.734 -0.406 -0.406 -0.734 0.281 0.557
## StPBc:TMHW2 0.534 -0.389 -0.389 -0.702 -0.406 -0.734 0.281 0.295
## SPB:TMHW1 SBH:TMHW2 SC:TMHW2
## SitBnnyHlls
## SiteCronull
## SitePalmBch
## TretmntMHW1
## TretmntMHW2
## StBHl:TMHW1
## StCrn:TMHW1
## StPBc:TMHW1
## StBHl:TMHW2 0.283
## StCrn:TMHW2 0.299 0.510
## StPBc:TMHW2 0.517 0.510 0.539
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')overdisp.glmer(my_glmer1)## Residual deviance: 1.255 on 102 degrees of freedom (ratio: 0.012)my_glmer_anova <- car::Anova(my_glmer1)
my_glmer_anova## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Fv_Fm
## Chisq Df Pr(>Chisq)
## Site 4.3046 3 0.23040
## Treatment 89.4045 2 < 2e-16 ***
## Site:Treatment 12.7934 6 0.04644 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- Post-hoc test
emm3 <- emmeans(my_glmer1, specs = pairwise ~ Treatment|Site, type = "response", adjust = "tukey")
summary(emm3)## $emmeans
## Site = Black Head:
## Treatment emmean SE df lower.CL upper.CL
## Ambient 0.2577 0.0394 98.8 0.1796 0.336
## MHW1 0.2785 0.0333 97.2 0.2123 0.345
## MHW2 0.0871 0.0370 97.4 0.0136 0.161
##
## Site = Bonny Hills:
## Treatment emmean SE df lower.CL upper.CL
## Ambient 0.3246 0.0369 99.5 0.2513 0.398
## MHW1 0.2605 0.0393 98.8 0.1824 0.339
## MHW2 0.0943 0.0420 100.3 0.0109 0.178
##
## Site = Cronulla:
## Treatment emmean SE df lower.CL upper.CL
## Ambient 0.3611 0.0370 97.4 0.2876 0.435
## MHW1 0.3784 0.0319 96.1 0.3152 0.442
## MHW2 0.0543 0.0333 97.2 -0.0118 0.120
##
## Site = Palm Beach:
## Treatment emmean SE df lower.CL upper.CL
## Ambient 0.2947 0.0333 97.2 0.2286 0.361
## MHW1 0.2916 0.0333 97.2 0.2255 0.358
## MHW2 0.1697 0.0369 99.5 0.0964 0.243
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## Site = Black Head:
## contrast estimate SE df t.ratio p.value
## Ambient - MHW1 -0.02070 0.0517 101.1 -0.401 0.9153
## Ambient - MHW2 0.17064 0.0539 100.7 3.165 0.0058
## MHW1 - MHW2 0.19134 0.0499 100.8 3.837 0.0006
##
## Site = Bonny Hills:
## contrast estimate SE df t.ratio p.value
## Ambient - MHW1 0.06406 0.0538 99.9 1.190 0.4619
## Ambient - MHW2 0.23027 0.0557 99.3 4.131 0.0002
## MHW1 - MHW2 0.16621 0.0575 101.0 2.891 0.0130
##
## Site = Cronulla:
## contrast estimate SE df t.ratio p.value
## Ambient - MHW1 -0.01731 0.0488 100.2 -0.354 0.9332
## Ambient - MHW2 0.30684 0.0497 99.2 6.175 <.0001
## MHW1 - MHW2 0.32414 0.0461 98.6 7.032 <.0001
##
## Site = Palm Beach:
## contrast estimate SE df t.ratio p.value
## Ambient - MHW1 0.00309 0.0471 99.1 0.066 0.9976
## Ambient - MHW2 0.12506 0.0497 99.0 2.517 0.0355
## MHW1 - MHW2 0.12197 0.0497 99.0 2.455 0.0415
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimatesemm3 <- emmeans(my_glmer1, specs = pairwise ~ Site|Treatment, type = "response", adjust = "tukey")
summary(emm3)## $emmeans
## Treatment = Ambient:
## Site emmean SE df lower.CL upper.CL
## Black Head 0.2577 0.0394 98.8 0.1796 0.336
## Bonny Hills 0.3246 0.0369 99.5 0.2513 0.398
## Cronulla 0.3611 0.0370 97.4 0.2876 0.435
## Palm Beach 0.2947 0.0333 97.2 0.2286 0.361
##
## Treatment = MHW1:
## Site emmean SE df lower.CL upper.CL
## Black Head 0.2785 0.0333 97.2 0.2123 0.345
## Bonny Hills 0.2605 0.0393 98.8 0.1824 0.339
## Cronulla 0.3784 0.0319 96.1 0.3152 0.442
## Palm Beach 0.2916 0.0333 97.2 0.2255 0.358
##
## Treatment = MHW2:
## Site emmean SE df lower.CL upper.CL
## Black Head 0.0871 0.0370 97.4 0.0136 0.161
## Bonny Hills 0.0943 0.0420 100.3 0.0109 0.178
## Cronulla 0.0543 0.0333 97.2 -0.0118 0.120
## Palm Beach 0.1697 0.0369 99.5 0.0964 0.243
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## Treatment = Ambient:
## contrast estimate SE df t.ratio p.value
## Black Head - Bonny Hills -0.06681 0.0537 98.8 -1.244 0.6005
## Black Head - Cronulla -0.10336 0.0538 99.8 -1.922 0.2256
## Black Head - Palm Beach -0.03698 0.0515 100.3 -0.718 0.8899
## Bonny Hills - Cronulla -0.03656 0.0521 99.3 -0.701 0.8963
## Bonny Hills - Palm Beach 0.02983 0.0497 99.0 0.600 0.9317
## Cronulla - Palm Beach 0.06638 0.0498 100.8 1.332 0.5450
##
## Treatment = MHW1:
## contrast estimate SE df t.ratio p.value
## Black Head - Bonny Hills 0.01795 0.0516 101.1 0.348 0.9855
## Black Head - Cronulla -0.09996 0.0461 98.6 -2.168 0.1395
## Black Head - Palm Beach -0.01318 0.0472 99.1 -0.280 0.9923
## Bonny Hills - Cronulla -0.11792 0.0506 100.4 -2.329 0.0983
## Bonny Hills - Palm Beach -0.03114 0.0516 101.1 -0.603 0.9308
## Cronulla - Palm Beach 0.08678 0.0461 98.6 1.882 0.2423
##
## Treatment = MHW2:
## contrast estimate SE df t.ratio p.value
## Black Head - Bonny Hills -0.00717 0.0557 99.0 -0.129 0.9992
## Black Head - Cronulla 0.03284 0.0497 99.2 0.661 0.9115
## Black Head - Palm Beach -0.08256 0.0521 99.3 -1.584 0.3925
## Bonny Hills - Cronulla 0.04001 0.0535 99.6 0.747 0.8776
## Bonny Hills - Palm Beach -0.07538 0.0557 99.3 -1.352 0.5322
## Cronulla - Palm Beach -0.11539 0.0497 99.0 -2.323 0.0998
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 4 estimates- Check model assumptions
res <- resid(my_glmer1)
#produce residual vs. fitted plot
plot(fitted(my_glmer1), res)
#add a horizontal line at 0
abline(0,0)
#create Q-Q plot for residuals
qqnorm(res)
#add a straight diagonal line to the plot
qqline(res) 
2.4. Day 22
- Run LMM for photosynthetic efficiency at day 22
my_glmer1<- lmer(Fv_Fm~ Site*Treatment + (1|Location), data=Juveniles_stat22_1)
summary(my_glmer1)## Linear mixed model fit by REML ['lmerMod']
## Formula: Fv_Fm ~ Site * Treatment + (1 | Location)
## Data: Juveniles_stat22_1
##
## REML criterion at convergence: -70.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.40305 -0.23558 0.07468 0.42078 2.50086
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 0.001963 0.04431
## Residual 0.015619 0.12497
## Number of obs: 79, groups: Location, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.2826151 0.0458021 6.170
## SiteBonny Hills 0.0503286 0.0590873 0.852
## SiteCronulla 0.0009258 0.0575091 0.016
## SitePalm Beach -0.0012438 0.0566122 -0.022
## TreatmentMHW1 -0.0139916 0.0581092 -0.241
## SiteBonny Hills:TreatmentMHW1 -0.1313075 0.0839806 -1.564
## SiteCronulla:TreatmentMHW1 0.0420341 0.0787639 0.534
## SitePalm Beach:TreatmentMHW1 -0.0316052 0.0798874 -0.396
##
## Correlation of Fixed Effects:
## (Intr) StBnnH StCrnl StPlmB TrMHW1 SBH:TM SC:TMH
## SitBnnyHlls -0.649
## SiteCronull -0.665 0.514
## SitePalmBch -0.685 0.527 0.539
## TretmntMHW1 -0.669 0.512 0.526 0.542
## StBHl:TMHW1 0.456 -0.699 -0.361 -0.367 -0.687
## StCrn:TMHW1 0.488 -0.375 -0.731 -0.395 -0.733 0.505
## StPBc:TMHW1 0.490 -0.376 -0.384 -0.713 -0.728 0.497 0.532overdisp.glmer(my_glmer1)## Residual deviance: 1.061 on 70 degrees of freedom (ratio: 0.015)my_glmer_anova <- car::Anova(my_glmer1)
my_glmer_anova## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Fv_Fm
## Chisq Df Pr(>Chisq)
## Site 1.4072 3 0.7038
## Treatment 1.9051 1 0.1675
## Site:Treatment 4.8201 3 0.1855- Check model assumptions
res <- resid(my_glmer1)
#produce residual vs. fitted plot
plot(fitted(my_glmer1), res)
#add a horizontal line at 0
abline(0,0)
#create Q-Q plot for residuals
qqnorm(res)
#add a straight diagonal line to the plot
qqline(res) 
3h) Length
1. Visualise differences
Length_means$Site <- factor(
Length_means$Site,
levels = c("Bonny Hills","Black Head","Palm Beach","Cronulla"),
labels = c(
"Bonny Hills","Black Head","Palm Beach","Cronulla"
))
length_plot <- Length_means %>%
filter(Treatment != "MHW2") %>%
ggplot(aes(x= Treatment, y=Mean_length, fill = Site)) +
geom_bar(stat="identity", color="black",
position=position_dodge()) +
geom_errorbar(aes(ymin=Mean_length-se, ymax=Mean_length+se), width=.2,
position=position_dodge(0.9))+
Theme1 +
labs( y=expression(paste("RGR(%d"^-1*")")), x="MHW scenario", fill = "Population") +
scale_fill_manual(values = c("darkred", "orange", "darkolivegreen3", "royalblue3")) +
scale_y_continuous(limits = c(-6, 3), breaks = seq(-7, 4, by = 2))
length_plot
ggsave("length_juv.jpeg", width = 140, height = 80, units = "mm")2. Statistical analyses
- Keep only alive individuals from MHW1 and ambient
Juveniles_length1 <- filter(Juveniles_length, Treatment != "MHW2")
Juveniles_length1 <- inner_join(Juveniles_length1, Juv_1, by=c("Code"))
Juveniles_length2 <- filter(Juveniles_length1, Survival != "DEAD")- Run LMM on relative growth rate (diff)
my_lmer <- lmer(sqrt(diff) ~ Site.y*Treatment.y + (1|Location), data=Juveniles_length2)## Warning in sqrt(diff): NaNs produced
## Warning in sqrt(diff): NaNs produced
## Warning in sqrt(diff): NaNs producedsummary(my_lmer)## Linear mixed model fit by REML ['lmerMod']
## Formula: sqrt(diff) ~ Site.y * Treatment.y + (1 | Location)
## Data: Juveniles_length2
##
## REML criterion at convergence: 15
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8349 -0.3606 -0.1070 0.5630 2.1325
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 0.01489 0.1220
## Residual 0.06794 0.2607
## Number of obs: 28, groups: Location, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.48745 0.14575 3.345
## Site.yBonny Hills 0.11028 0.23705 0.465
## Site.yCronulla -0.03022 0.17791 -0.170
## Site.yPalm Beach 0.21452 0.23687 0.906
## Treatment.yMHW1 0.14323 0.23994 0.597
## Site.yBonny Hills:Treatment.yMHW1 -0.37396 0.33236 -1.125
## Site.yCronulla:Treatment.yMHW1 0.05631 0.29087 0.194
## Site.yPalm Beach:Treatment.yMHW1 -0.33551 0.35240 -0.952
##
## Correlation of Fixed Effects:
## (Intr) St.yBH St.yCr St.yPB T.MHW1 S.BH:T S.C:T.
## St.yBnnyHll -0.553
## Site.yCrnll -0.691 0.439
## Sit.yPlmBch -0.547 0.358 0.421
## Trtmnt.MHW1 -0.558 0.353 0.432 0.315
## S.BH:T.MHW1 0.395 -0.692 -0.310 -0.223 -0.728
## St.C:T.MHW1 0.452 -0.286 -0.623 -0.259 -0.831 0.612
## S.PB:T.MHW1 0.388 -0.269 -0.287 -0.681 -0.685 0.510 0.583my_lmer_anova <- car::Anova(my_lmer)
my_lmer_anova## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: sqrt(diff)
## Chisq Df Pr(>Chisq)
## Site.y 1.2400 3 0.7434
## Treatment.y 0.1030 1 0.7483
## Site.y:Treatment.y 3.4395 3 0.3287- Check model assumptions
res <- resid(my_lmer)
#produce residual vs. fitted plot
plot(fitted(my_lmer), res)
#add a horizontal line at 0
abline(0,0)
#create Q-Q plot for residuals
qqnorm(res)
#add a straight diagonal line to the plot
qqline(res) 
3i) C:N ratio
1. Visualise differences
CNresults_mean$Site.x <- factor(
CNresults_mean$Site.x,
levels = c("Bonny Hills","Black Head","Palm Beach","Cronulla"),
labels = c(
"Bonny Hills","Black Head","Palm Beach","Cronulla"
))ggplot(data =CNresults_mean, aes(x = Treatment.y, y =Mean_CNratio, fill = Site.x)) +
geom_bar(stat="identity", color="black",
position=position_dodge()) +
geom_errorbar(aes(ymin=Mean_CNratio-se, ymax=Mean_CNratio+se), width=.2,
position=position_dodge(0.9))+
scale_fill_manual(values = c("darkred","orange", "darkolivegreen3", "royalblue3")) +
labs( y=expression(paste("C:N ratio")), x="MHW scenario", fill = "Population") +
Theme1
ggsave("CN.jpeg", width = 150, height = 70, units = "mm")2. Statistical analyses
- Run LMM for CN ratio
CN_results5 <- inner_join(CN_results4 , Tanks_metadata, by=c("Code"))
my_lmer <- lmer(sqrt(CN_ratio) ~ Site.x*Treatment.y + (1|Location), data=CN_results5)## boundary (singular) fit: see help('isSingular')summary(my_lmer)## Linear mixed model fit by REML ['lmerMod']
## Formula: sqrt(CN_ratio) ~ Site.x * Treatment.y + (1 | Location)
## Data: CN_results5
##
## REML criterion at convergence: 80.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.7445 -0.6990 0.1454 0.6770 1.6718
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 0.000 0.0000
## Residual 0.339 0.5822
## Number of obs: 46, groups: Location, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 6.1246 0.2377 25.768
## Site.xBonny Hills -0.2310 0.3361 -0.687
## Site.xCronulla -0.1626 0.3361 -0.484
## Site.xPalm Beach 0.4657 0.3361 1.385
## Treatment.yMHW1 0.5020 0.3361 1.493
## Site.xBonny Hills:Treatment.yMHW1 0.2770 0.4871 0.569
## Site.xCronulla:Treatment.yMHW1 -0.5748 0.4754 -1.209
## Site.xPalm Beach:Treatment.yMHW1 -0.8416 0.4871 -1.728
##
## Correlation of Fixed Effects:
## (Intr) St.xBH St.xCr St.xPB T.MHW1 S.BH:T S.C:T.
## St.xBnnyHll -0.707
## Site.xCrnll -0.707 0.500
## Sit.xPlmBch -0.707 0.500 0.500
## Trtmnt.MHW1 -0.707 0.500 0.500 0.500
## S.BH:T.MHW1 0.488 -0.690 -0.345 -0.345 -0.690
## St.C:T.MHW1 0.500 -0.354 -0.707 -0.354 -0.707 0.488
## S.PB:T.MHW1 0.488 -0.345 -0.345 -0.690 -0.690 0.476 0.488
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')my_lmer_anova <- car::Anova(my_lmer)
my_lmer_anova## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: sqrt(CN_ratio)
## Chisq Df Pr(>Chisq)
## Site.x 5.5239 3 0.13722
## Treatment.y 1.5911 1 0.20717
## Site.x:Treatment.y 6.4958 3 0.08983 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- Post-hoc test
emm3 = emmeans(my_lmer, specs = pairwise ~ Site.x, type = "response", adjust = "none")## NOTE: Results may be misleading due to involvement in interactionsemm3## $emmeans
## Site.x response SE df lower.CL upper.CL
## Black Head 40.6 2.22 25.6 36.2 45.3
## Bonny Hills 39.5 2.28 24.5 34.9 44.3
## Cronulla 35.1 2.03 28.8 31.1 39.4
## Palm Beach 41.2 2.31 28.2 36.6 46.1
##
## Results are averaged over the levels of: Treatment.y
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Intervals are back-transformed from the sqrt scale
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## Black Head - Bonny Hills 0.0925 0.247 36.4 0.375 0.7100
## Black Head - Cronulla 0.4499 0.240 35.3 1.873 0.0693
## Black Head - Palm Beach -0.0449 0.247 35.8 -0.181 0.8571
## Bonny Hills - Cronulla 0.3574 0.248 37.2 1.440 0.1582
## Bonny Hills - Palm Beach -0.1374 0.252 35.4 -0.546 0.5886
## Cronulla - Palm Beach -0.4948 0.246 36.0 -2.011 0.0519
##
## Results are averaged over the levels of: Treatment.y
## Note: contrasts are still on the sqrt scale
## Degrees-of-freedom method: kenward-roger- Check model assumptions
res <- resid(my_lmer)
#produce residual vs. fitted plot
plot(fitted(my_lmer), res)
#add a horizontal line at 0
abline(0,0)
#create Q-Q plot for residuals
qqnorm(res)
#add a straight diagonal line to the plot
qqline(res) 
3j) Carbon content
1. Visualise differences
ggplot(data =CNresults_mean, aes(x = Treatment.y, y =Mean_C, fill = Site.x)) +
geom_bar(stat="identity", color="black",
position=position_dodge()) +
geom_errorbar(aes(ymin=Mean_C-se_C, ymax=Mean_C+se_C), width=.2,
position=position_dodge(0.9))+
scale_fill_manual(values = c("darkred","orange", "darkolivegreen3", "royalblue3")) +
labs( y=expression(paste("Tissue Carbon (%)")), x="MHW scenario", fill = "Population") +
Theme1
ggsave("Carbon.jpeg", width = 150, height = 70, units = "mm")2. Statistical analyses
- Run LMM for carbon content
my_lmer <- lmer(sqrt(C) ~ Site.x*Treatment.y + (1|Location), data=CN_results5)## boundary (singular) fit: see help('isSingular')summary(my_lmer)## Linear mixed model fit by REML ['lmerMod']
## Formula: sqrt(C) ~ Site.x * Treatment.y + (1 | Location)
## Data: CN_results5
##
## REML criterion at convergence: -17.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.53125 -0.56353 -0.02449 0.62616 2.04468
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 0.00000 0.0000
## Residual 0.02539 0.1593
## Number of obs: 46, groups: Location, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 5.60043 0.06505 86.097
## Site.xBonny Hills 0.03056 0.09199 0.332
## Site.xCronulla 0.16494 0.09199 1.793
## Site.xPalm Beach 0.27927 0.09199 3.036
## Treatment.yMHW1 0.08698 0.09199 0.945
## Site.xBonny Hills:Treatment.yMHW1 -0.03534 0.13331 -0.265
## Site.xCronulla:Treatment.yMHW1 -0.10017 0.13010 -0.770
## Site.xPalm Beach:Treatment.yMHW1 -0.11961 0.13331 -0.897
##
## Correlation of Fixed Effects:
## (Intr) St.xBH St.xCr St.xPB T.MHW1 S.BH:T S.C:T.
## St.xBnnyHll -0.707
## Site.xCrnll -0.707 0.500
## Sit.xPlmBch -0.707 0.500 0.500
## Trtmnt.MHW1 -0.707 0.500 0.500 0.500
## S.BH:T.MHW1 0.488 -0.690 -0.345 -0.345 -0.690
## St.C:T.MHW1 0.500 -0.354 -0.707 -0.354 -0.707 0.488
## S.PB:T.MHW1 0.488 -0.345 -0.345 -0.690 -0.690 0.476 0.488
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')my_lmer_anova <- car::Anova(my_lmer)
my_lmer_anova## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: sqrt(C)
## Chisq Df Pr(>Chisq)
## Site.x 14.3015 3 0.002522 **
## Treatment.y 0.2566 1 0.612462
## Site.x:Treatment.y 1.0587 3 0.787058
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- Post-hoc test
emm3 = emmeans(my_lmer, specs = pairwise ~ Site.x, type = "response", adjust = "none")## NOTE: Results may be misleading due to involvement in interactionsemm3## $emmeans
## Site.x response SE df lower.CL upper.CL
## Black Head 31.9 0.538 25.6 30.8 33.0
## Bonny Hills 32.0 0.562 24.5 30.9 33.2
## Cronulla 33.2 0.539 28.8 32.1 34.3
## Palm Beach 34.4 0.577 28.2 33.2 35.6
##
## Results are averaged over the levels of: Treatment.y
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Intervals are back-transformed from the sqrt scale
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## Black Head - Bonny Hills -0.0129 0.0676 36.4 -0.191 0.8498
## Black Head - Cronulla -0.1149 0.0657 35.3 -1.747 0.0893
## Black Head - Palm Beach -0.2195 0.0677 35.8 -3.243 0.0026
## Bonny Hills - Cronulla -0.1020 0.0679 37.2 -1.501 0.1418
## Bonny Hills - Palm Beach -0.2066 0.0689 35.4 -3.000 0.0049
## Cronulla - Palm Beach -0.1046 0.0673 36.0 -1.553 0.1291
##
## Results are averaged over the levels of: Treatment.y
## Note: contrasts are still on the sqrt scale
## Degrees-of-freedom method: kenward-roger- Check model assumptions
res <- resid(my_lmer)
#produce residual vs. fitted plot
plot(fitted(my_lmer), res)
#add a horizontal line at 0
abline(0,0)
#create Q-Q plot for residuals
qqnorm(res)
#add a straight diagonal line to the plot
qqline(res) 
3k) Nitrogen content
1. Visualise differences
ggplot(data =CNresults_mean, aes(x = Treatment.y, y =Mean_N, fill = Site.x)) +
geom_bar(stat="identity", color="black",
position=position_dodge()) +
geom_errorbar(aes(ymin=Mean_N-se_N, ymax=Mean_N+se_N), width=.2,
position=position_dodge(0.9))+
scale_fill_manual(values = c("darkred","orange", "darkolivegreen3", "royalblue3")) +
labs( y=expression(paste("Tissue Nitrogen (%)")), x="MHW scenario", fill = "Population") +
Theme1
ggsave("Nitrogen.jpeg", width = 150, height = 70, units = "mm")2. Statistical analyses
- Run LMM for nitrogen content
my_lmer <- lmer(sqrt(N) ~ Treatment.y*Site.x + (1|Location), data=CN_results5)## boundary (singular) fit: see help('isSingular')summary(my_lmer)## Linear mixed model fit by REML ['lmerMod']
## Formula: sqrt(N) ~ Treatment.y * Site.x + (1 | Location)
## Data: CN_results5
##
## REML criterion at convergence: -74.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.3832 -0.8255 -0.1673 0.5817 1.6883
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 0.000000 0.00000
## Residual 0.005702 0.07551
## Number of obs: 46, groups: Location, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.91838 0.03083 29.792
## Treatment.yMHW1 -0.05087 0.04360 -1.167
## Site.xBonny Hills 0.04012 0.04360 0.920
## Site.xCronulla 0.05320 0.04360 1.220
## Site.xPalm Beach -0.01891 0.04360 -0.434
## Treatment.yMHW1:Site.xBonny Hills -0.04819 0.06318 -0.763
## Treatment.yMHW1:Site.xCronulla 0.05836 0.06165 0.947
## Treatment.yMHW1:Site.xPalm Beach 0.09028 0.06318 1.429
##
## Correlation of Fixed Effects:
## (Intr) Tr.MHW1 St.xBH St.xCr St.xPB T.MHWH T.MHW1:
## Trtmnt.MHW1 -0.707
## St.xBnnyHll -0.707 0.500
## Site.xCrnll -0.707 0.500 0.500
## Sit.xPlmBch -0.707 0.500 0.500 0.500
## T.MHW1:S.BH 0.488 -0.690 -0.690 -0.345 -0.345
## Tr.MHW1:S.C 0.500 -0.707 -0.354 -0.707 -0.354 0.488
## T.MHW1:S.PB 0.488 -0.690 -0.345 -0.345 -0.690 0.476 0.488
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')my_lmer_anova <- car::Anova(my_lmer)
my_lmer_anova## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: sqrt(N)
## Chisq Df Pr(>Chisq)
## Treatment.y 1.3135 1 0.25176
## Site.x 7.9505 3 0.04705 *
## Treatment.y:Site.x 5.5152 3 0.13773
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- Post-hoc test
emm3 = emmeans(my_lmer, specs = pairwise ~ Site.x, type = "response", adjust = "none")## NOTE: Results may be misleading due to involvement in interactionsemm3## $emmeans
## Site.x response SE df lower.CL upper.CL
## Black Head 0.797 0.0404 25.6 0.716 0.883
## Bonny Hills 0.826 0.0428 24.5 0.740 0.917
## Cronulla 0.951 0.0433 28.8 0.865 1.042
## Palm Beach 0.845 0.0429 28.2 0.759 0.935
##
## Results are averaged over the levels of: Treatment.y
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Intervals are back-transformed from the sqrt scale
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## Black Head - Bonny Hills -0.0160 0.0320 36.4 -0.500 0.6198
## Black Head - Cronulla -0.0824 0.0312 35.3 -2.645 0.0121
## Black Head - Palm Beach -0.0262 0.0321 35.8 -0.818 0.4188
## Bonny Hills - Cronulla -0.0664 0.0322 37.2 -2.061 0.0463
## Bonny Hills - Palm Beach -0.0102 0.0326 35.4 -0.313 0.7562
## Cronulla - Palm Beach 0.0562 0.0319 36.0 1.759 0.0870
##
## Results are averaged over the levels of: Treatment.y
## Note: contrasts are still on the sqrt scale
## Degrees-of-freedom method: kenward-roger- Check model assumptions
res <- resid(my_lmer)
#produce residual vs. fitted plot
plot(fitted(my_lmer), res)
#add a horizontal line at 0
abline(0,0)
#create Q-Q plot for residuals
qqnorm(res)
#add a straight diagonal line to the plot
qqline(res) 
3l) Chlorophyll a
1. Visualise differences
Pigments_mean$Site_Code.x <- factor(
Pigments_mean$Site_Code.x,
levels = c("Bonny Hills","Black Head","Palm Beach","Cronulla"),
labels = c("Bonny Hills","Black Head","Palm Beach","Cronulla"
))
Pigments_mean4 <- filter(Pigments_mean, Treatment.y != "Initial")ggplot(data =Pigments_mean4, aes(x = Treatment.y, y =Mean_chla, fill = Site_Code.x)) +
geom_bar(stat="identity", color="black",
position=position_dodge()) +
geom_errorbar(aes(ymin=Mean_chla-se_chla, ymax=Mean_chla+se_chla), width=.2,
position=position_dodge(0.9))+
scale_fill_manual(values = c("darkred","orange", "darkolivegreen3", "royalblue3")) +
labs( y=expression(paste("Chlorophyll a (mg g"^-1*")")), x="MHW scenario", fill = "Population") +
theme_classic()
ggsave("chla.jpeg", width = 150, height = 70, units = "mm")2. Statistical analyses
- Run LMM for chlorophyll-a
Pigments7 <- inner_join(Pigments6 , Tanks_metadata, by=c("Code"))
my_lmer <- lmer(sqrt(chla) ~ Site_Code.x*Treatment.y + (1|Location), data=Pigments7) # Location is shelf## boundary (singular) fit: see help('isSingular')summary(my_lmer)## Linear mixed model fit by REML ['lmerMod']
## Formula: sqrt(chla) ~ Site_Code.x * Treatment.y + (1 | Location)
## Data: Pigments7
##
## REML criterion at convergence: -216.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.45970 -0.69387 -0.07341 0.60190 2.59818
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 0.000000 0.00000
## Residual 0.000421 0.02052
## Number of obs: 55, groups: Location, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.089480 0.008377 10.682
## Site_Code.xBonny Hills -0.001068 0.011081 -0.096
## Site_Code.xCronulla 0.001659 0.011416 0.145
## Site_Code.xPalm Beach 0.002309 0.010814 0.214
## Treatment.yMHW1 0.014446 0.011847 1.219
## Site_Code.xBonny Hills:Treatment.yMHW1 -0.008125 0.016648 -0.488
## Site_Code.xCronulla:Treatment.yMHW1 -0.005035 0.015910 -0.317
## Site_Code.xPalm Beach:Treatment.yMHW1 -0.015115 0.016040 -0.942
##
## Correlation of Fixed Effects:
## (Intr) St_C.BH St_C.C St_C.PB T.MHW1 S_C.BH: S_C.C:
## St_Cd.xBnnH -0.756
## St_Cd.xCrnl -0.734 0.555
## St_Cd.xPlmB -0.775 0.586 0.568
## Trtmnt.MHW1 -0.707 0.535 0.519 0.548
## S_C.BH:T.MH 0.503 -0.666 -0.369 -0.390 -0.712
## S_C.C:T.MHW 0.527 -0.398 -0.718 -0.408 -0.745 0.530
## S_C.PB:T.MH 0.522 -0.395 -0.383 -0.674 -0.739 0.526 0.550
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')my_lmer_anova <- car::Anova(my_lmer)
my_lmer_anova## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: sqrt(chla)
## Chisq Df Pr(>Chisq)
## Site_Code.x 0.6091 3 0.8943
## Treatment.y 1.6159 1 0.2037
## Site_Code.x:Treatment.y 0.9524 3 0.8128- Check model assumptions
res <- resid(my_lmer)
#produce residual vs. fitted plot
plot(fitted(my_lmer), res)
#add a horizontal line at 0
abline(0,0)
#create Q-Q plot for residuals
qqnorm(res)
#add a straight diagonal line to the plot
qqline(res) 
3m) Chlorophyll c
1. Visualise differences
ggplot(data =Pigments_mean, aes(x = Treatment.y, y =Mean_chlc, fill = Site_Code.x)) +
geom_bar(stat="identity", color="black",
position=position_dodge()) +
geom_errorbar(aes(ymin=Mean_chlc-se_chlc, ymax=Mean_chlc+se_chlc), width=.2,
position=position_dodge(0.9))+
scale_fill_manual(values = c("darkred","orange", "darkolivegreen3", "royalblue3")) +
labs( y=expression(paste("Chlorophyll c (mg g"^-1*")")), x="MHW scenario", fill = "Population") +
Theme1
ggsave("chlc.jpeg", width = 150, height = 70, units = "mm")2. Statistical analyses
- Run LMM for chlorophyll-c
my_lmer <- lmer(sqrt(chlc) ~ Site_Code.x*Treatment.y + (1|Location), data=Pigments7)
summary(my_lmer)## Linear mixed model fit by REML ['lmerMod']
## Formula: sqrt(chlc) ~ Site_Code.x * Treatment.y + (1 | Location)
## Data: Pigments7
##
## REML criterion at convergence: -365.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.0758 -0.5949 0.0172 0.5566 2.3716
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 9.915e-07 0.0009958
## Residual 1.697e-05 0.0041198
## Number of obs: 55, groups: Location, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.0569419 0.0017531 32.480
## Site_Code.xBonny Hills -0.0003781 0.0022305 -0.170
## Site_Code.xCronulla -0.0007498 0.0022949 -0.327
## Site_Code.xPalm Beach -0.0008956 0.0021836 -0.410
## Treatment.yMHW1 0.0037292 0.0024006 1.553
## Site_Code.xBonny Hills:Treatment.yMHW1 -0.0002664 0.0033573 -0.079
## Site_Code.xCronulla:Treatment.yMHW1 0.0022075 0.0032029 0.689
## Site_Code.xPalm Beach:Treatment.yMHW1 -0.0019781 0.0032441 -0.610
##
## Correlation of Fixed Effects:
## (Intr) St_C.BH St_C.C St_C.PB T.MHW1 S_C.BH: S_C.C:
## St_Cd.xBnnH -0.731
## St_Cd.xCrnl -0.705 0.553
## St_Cd.xPlmB -0.752 0.588 0.565
## Trtmnt.MHW1 -0.685 0.536 0.515 0.550
## S_C.BH:T.MH 0.483 -0.662 -0.369 -0.386 -0.710
## S_C.C:T.MHW 0.504 -0.397 -0.717 -0.406 -0.742 0.531
## S_C.PB:T.MH 0.504 -0.397 -0.382 -0.674 -0.740 0.526 0.552my_lmer_anova <- car::Anova(my_lmer)
my_lmer_anova## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: sqrt(chlc)
## Chisq Df Pr(>Chisq)
## Site_Code.x 2.1338 3 0.54511
## Treatment.y 10.7194 1 0.00106 **
## Site_Code.x:Treatment.y 1.9035 3 0.59268
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- Post-hoc test
emm3 = emmeans(my_lmer, specs = pairwise ~ Treatment.y, type = "response", adjust = "none")## NOTE: Results may be misleading due to involvement in interactionsemm3## $emmeans
## Treatment.y response SE df lower.CL upper.CL
## Ambient 0.00319 0.000103 8.32 0.00295 0.00343
## MHW1 0.00362 0.000116 12.14 0.00337 0.00388
##
## Results are averaged over the levels of: Site_Code.x
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Intervals are back-transformed from the sqrt scale
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## Ambient - MHW1 -0.00372 0.00117 46.4 -3.170 0.0027
##
## Results are averaged over the levels of: Site_Code.x
## Note: contrasts are still on the sqrt scale
## Degrees-of-freedom method: kenward-roger- Check model assumptions
res <- resid(my_lmer)
#produce residual vs. fitted plot
plot(fitted(my_lmer), res)
#add a horizontal line at 0
abline(0,0)
#create Q-Q plot for residuals
qqnorm(res)
#add a straight diagonal line to the plot
qqline(res) 
3n) Fucoxanthin
1. Visualise differences
ggplot(data =Pigments_mean, aes(x = Treatment.y, y =Mean_Fucoxanthin, fill = Site_Code.x)) +
geom_bar(stat="identity", color="black",
position=position_dodge()) +
geom_errorbar(aes(ymin=Mean_Fucoxanthin-se_Fucoxanthin, ymax=Mean_Fucoxanthin+se_Fucoxanthin), width=.2,
position=position_dodge(0.9))+
scale_fill_manual(values = c("darkred","orange", "darkolivegreen3", "royalblue3")) +
labs( y=expression(paste("Fucoxanthin (mg g"^-1*")")), x="MHW scenario", fill = "Population") +
Theme1
ggsave("fucoxanthin.jpeg", width = 150, height = 70, units = "mm")2. Statistical analyses
- Run LMM for fucoxanthins
my_lmer <- lmer(sqrt(Fucoxanthin) ~ Site_Code.x*Treatment.y + (1|Location), data=Pigments7)## boundary (singular) fit: see help('isSingular')summary(my_lmer)## Linear mixed model fit by REML ['lmerMod']
## Formula: sqrt(Fucoxanthin) ~ Site_Code.x * Treatment.y + (1 | Location)
## Data: Pigments7
##
## REML criterion at convergence: -229.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.0289 -0.6197 -0.1589 0.5990 2.2257
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 0.0000000 0.00000
## Residual 0.0003184 0.01784
## Number of obs: 55, groups: Location, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 5.077e-02 7.284e-03 6.970
## Site_Code.xBonny Hills -4.167e-03 9.636e-03 -0.432
## Site_Code.xCronulla -7.651e-04 9.927e-03 -0.077
## Site_Code.xPalm Beach 7.809e-05 9.404e-03 0.008
## Treatment.yMHW1 1.194e-02 1.030e-02 1.159
## Site_Code.xBonny Hills:Treatment.yMHW1 -7.669e-03 1.448e-02 -0.530
## Site_Code.xCronulla:Treatment.yMHW1 -7.634e-03 1.383e-02 -0.552
## Site_Code.xPalm Beach:Treatment.yMHW1 -1.335e-02 1.395e-02 -0.957
##
## Correlation of Fixed Effects:
## (Intr) St_C.BH St_C.C St_C.PB T.MHW1 S_C.BH: S_C.C:
## St_Cd.xBnnH -0.756
## St_Cd.xCrnl -0.734 0.555
## St_Cd.xPlmB -0.775 0.586 0.568
## Trtmnt.MHW1 -0.707 0.535 0.519 0.548
## S_C.BH:T.MH 0.503 -0.666 -0.369 -0.390 -0.712
## S_C.C:T.MHW 0.527 -0.398 -0.718 -0.408 -0.745 0.530
## S_C.PB:T.MH 0.522 -0.395 -0.383 -0.674 -0.739 0.526 0.550
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')my_lmer_anova <- car::Anova(my_lmer)
my_lmer_anova## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: sqrt(Fucoxanthin)
## Chisq Df Pr(>Chisq)
## Site_Code.x 1.3414 3 0.7193
## Treatment.y 0.8428 1 0.3586
## Site_Code.x:Treatment.y 0.9173 3 0.8212- Check model assumptions
res <- resid(my_lmer)
#produce residual vs. fitted plot
plot(fitted(my_lmer), res)
#add a horizontal line at 0
abline(0,0)
#create Q-Q plot for residuals
qqnorm(res)
#add a straight diagonal line to the plot
qqline(res) 
R Session Information
In this part, information about the operating system and R packages attached during the production of this document can be found
sessionInfo() %>% panderR version 4.4.1 (2024-06-14 ucrt)
Platform: x86_64-w64-mingw32/x64
locale: LC_COLLATE=English_Australia.utf8, LC_CTYPE=English_Australia.utf8, LC_MONETARY=English_Australia.utf8, LC_NUMERIC=C and LC_TIME=English_Australia.utf8
attached base packages: stats, graphics, grDevices, utils, datasets, methods and base
other attached packages: pander(v.0.6.5), RVAideMemoire(v.0.9-83-7), gt(v.1.0.0), skimr(v.2.1.5), outliers(v.0.15), car(v.3.1-2), carData(v.3.0-5), emmeans(v.1.10.2), lme4(v.1.1-35.4), Matrix(v.1.7-0), ggplot2(v.3.5.1), patchwork(v.1.3.0), dplyr(v.1.1.4), tidyr(v.1.3.1) and summarytools(v.1.0.1)
loaded via a namespace (and not attached): tidyselect(v.1.2.1), farver(v.2.1.2), fastmap(v.1.2.0), TH.data(v.1.1-2), digest(v.0.6.35), estimability(v.1.5.1), timechange(v.0.3.0), lifecycle(v.1.0.4), survival(v.3.6-4), magrittr(v.2.0.3), compiler(v.4.4.1), rlang(v.1.1.4), sass(v.0.4.9), tools(v.4.4.1), yaml(v.2.3.8), knitr(v.1.50), labeling(v.0.4.3), plyr(v.1.8.9), xml2(v.1.3.6), repr(v.1.1.7), multcomp(v.1.4-26), abind(v.1.4-5), withr(v.3.0.2), purrr(v.1.0.2), grid(v.4.4.1), xtable(v.1.8-4), colorspace(v.2.1-0), scales(v.1.3.0), MASS(v.7.3-60.2), cli(v.3.6.2), mvtnorm(v.1.2-5), rmarkdown(v.2.29), ragg(v.1.3.2), generics(v.0.1.3), rstudioapi(v.0.16.0), reshape2(v.1.4.4), minqa(v.1.2.7), cachem(v.1.1.0), stringr(v.1.5.1), splines(v.4.4.1), parallel(v.4.4.1), matrixStats(v.1.5.0), base64enc(v.0.1-3), vctrs(v.0.6.5), boot(v.1.3-30), sandwich(v.3.1-1), jsonlite(v.1.9.1), rapportools(v.1.1), pbkrtest(v.0.5.2), systemfonts(v.1.2.1), magick(v.2.8.3), jquerylib(v.0.1.4), glue(v.1.7.0), nloptr(v.2.1.0), codetools(v.0.2-20), lubridate(v.1.9.3), stringi(v.1.8.4), gtable(v.0.3.6), munsell(v.0.5.1), tibble(v.3.2.1), pillar(v.1.10.1), htmltools(v.0.5.8.1), R6(v.2.6.1), tcltk(v.4.4.1), textshaping(v.0.4.0), evaluate(v.1.0.3), lattice(v.0.22-6), backports(v.1.5.0), broom(v.1.0.6), pryr(v.0.1.6), bslib(v.0.9.0), Rcpp(v.1.0.12), coda(v.0.19-4.1), nlme(v.3.1-167), checkmate(v.2.3.1), xfun(v.0.51), zoo(v.1.8-12) and pkgconfig(v.2.0.3)