Juvenile~AE 2026 - Data analysis
Water chemistry
QA/QC
Each sample is run in replicate on the Metrhom titrator. If the replicates are within an acceptable threshold (usually < 5umol difference), then the mean of the replicates is used as the final TA value for that sample. If not, the sample is re-run. If the delta is still >5 umol, sample will be removed from dataset.
Frequency of replicate deltas. All samples below 5 umol are acceptable
Sample replicate deltas broken down by treatment.
Remove all samples with delta >5 umol. Print eliminated samples.
Code
# Remove bottle samples with >5 umol delta
bottle_df <- bottle_j |>
filter(ta_diff<= 5 | is.na(ta_diff),
date >= "2026-01-01")
# Print which samples were removed
dplyr::anti_join(bottle_j, bottle_df, by = "bottle_name")# A tibble: 8 × 29
date parameter bottle_name tank treatment time salinity
<date> <chr> <chr> <fct> <fct> <dttm> <dbl>
1 2025-12-17 TA 14334 B1 ambient 1899-12-31 02:34:00 NA
2 2025-12-17 TA 14336 B3 high 1899-12-31 02:52:00 NA
3 2025-12-17 TA 14335 B2 medium 1899-12-31 02:42:00 NA
4 2025-12-18 TA 14337 B1 ambient 1899-12-31 04:07:00 NA
5 2025-12-18 TA 14339 B3 high 1899-12-31 04:19:00 NA
6 2025-12-18 TA 14338 B2 medium 1899-12-31 04:15:00 NA
7 2026-01-07 TA 14364 D3 high NA NA
8 2026-02-16 TA 14090 D2 medium 1899-12-31 10:00:00 35.2
# ℹ 22 more variables: temp <dbl>, glass_p_h <dbl>, dkh <dbl>, hanna_alk <dbl>,
# sw_dose <dbl>, ae_dose <dbl>, ro_dose <dbl>, interval <dbl>,
# diff_btwn_treatments <dbl>, notes <chr>, final_ta <dbl>, ta_diff <dbl>,
# batch <chr>, ta1_timestamp <dttm>, ta1 <dbl>, ta2_timestamp <dttm>,
# ta2 <dbl>, extra1_timestamp <dttm>, extra1_ta <dbl>,
# extra2_timestamp <dttm>, extra2_ta <dbl>, dic <dbl>Instrument precision
Get Metrohm robo TA instrument precision, i.e. SD of CRMs over all runs
Code
# Extract CRM runs from all rounds
# Define path to excel file with all roboTA rounds on different sheets
file_path <- "data/raw/metrohmOutput-2.xlsx"
# Retreive character vector of all sheet names in file
sheet_names <- excel_sheets(file_path)
# Create CRM dataframe
crm_all <- purrr::map_dfr(sheet_names[str_detect(sheet_names, "BAE")], ~ {
read_excel(file_path, sheet = .x) %>%
dplyr::filter(stringr::str_starts(`Bottle Name`, "w CRM 219")) %>%
mutate(
sheet = .x,
Timestamp = as.POSIXct(Timestamp, tz = "UTC") # force consistent type
)
}, .id = NULL) %>%
janitor::clean_names()[1] "SD = ± 2.84 umol/kg"[1] "n = 16"
Treatments
stats
Treatment Summary
mixed linear model.
Run mixed linear model to check for treatment significance and tank effects:
treat_lmer <- lmer(final_ta ~ treatment + (1|tank), data = bottle_df)
summary(treat_lmer)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: final_ta ~ treatment + (1 | tank)
Data: bottle_df
REML criterion at convergence: 1495
Scaled residuals:
Min 1Q Median 3Q Max
-1.7344 -0.3363 -0.0454 0.2085 9.1780
Random effects:
Groups Name Variance Std.Dev.
tank (Intercept) 1541 39.26
Residual 16495 128.43
Number of obs: 121, groups: tank, 6
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 2427.782 34.252 2.945 70.881 7.45e-06 ***
treatmentmedium 363.295 48.541 2.970 7.484 0.00512 **
treatmenthigh 781.937 48.548 2.971 16.106 0.00055 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) trtmntm
treatmntmdm -0.706
treatmnthgh -0.706 0.498 Random effects:
tank accounts for relatively small variance compared to residual.. the intraclass correlation (proportion of variance due to tank) is 1541/ (1541 + 16495) = 0.08.
So 8% of variance is attributable to tank differences.
Fixed effects: Treatment highly significant, but degrees of freedom is 2.9 kind of low.
Residuals: Max: 9.178. This is large (values beyond 3 are worth investigating). Points to an outlier.
Search for outlier:
Code
# Residuals vs fitted
plot(treat_lmer)
Code
# QQ plot
qqnorm(resid(treat_lmer))
qqline(resid(treat_lmer))
Likely in tank D3 during the period of doser malfunction. Lead to drop in salinity and spike in TA. Remove by filtering out tank samples from known malfunction period.
Code
bottle_join <- bottle_df |>
filter(!(tank == "D3" & date >="2026-02-04" & date <= "2026-02-05"))Re-run stats with outlier removed
treat_lmer <- lmer(final_ta ~ treatment + (1|tank), data = bottle_join)
summary(treat_lmer)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: final_ta ~ treatment + (1 | tank)
Data: bottle_join
REML criterion at convergence: 1323.1
Scaled residuals:
Min 1Q Median 3Q Max
-3.2748 -0.4501 -0.0101 0.5261 5.8266
Random effects:
Groups Name Variance Std.Dev.
tank (Intercept) 319.2 17.87
Residual 4240.7 65.12
Number of obs: 120, groups: tank, 6
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 2427.787 16.219 2.879 149.69 1.08e-06 ***
treatmentmedium 363.289 22.993 2.907 15.80 0.000656 ***
treatmenthigh 749.979 23.060 2.939 32.52 7.49e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) trtmntm
treatmntmdm -0.705
treatmnthgh -0.703 0.496 # Residuals vs fitted
plot(treat_lmer)
# QQ plot
qqnorm(resid(treat_lmer))
qqline(resid(treat_lmer))
tank variance down by 80%, residual variance down by 70%, max scaled residuals down to 5.6. Solid, move forward with other stats
ANOVA
Code
treat_aov <- aov(final_ta ~ treatment, data = bottle_join)
summary(treat_aov) Df Sum Sq Mean Sq F value Pr(>F)
treatment 2 11203939 5601969 1273 <2e-16 ***
Residuals 117 514993 4402
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
58 observations deleted due to missingnessCode
# Post-hoc comparison
treat_tukey <- glht(treat_aov,
linfct = mcp(treatment = "Tukey"))
summary(treat_tukey)
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Tukey Contrasts
Fit: aov(formula = final_ta ~ treatment, data = bottle_join)
Linear Hypotheses:
Estimate Std. Error t value Pr(>|t|)
medium - ambient == 0 363.22 14.74 24.64 <2e-16 ***
high - ambient == 0 748.69 14.84 50.45 <2e-16 ***
high - medium == 0 385.47 14.93 25.82 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)Treatments are sig different based on ANOVA.
Alkalinity (umol/kg) treatment table:
treatment final_ta.V1 final_ta.sd
1 ambient 2427.85 30.11
2 medium 2791.08 69.81
3 high 3176.55 87.16Code
# Quick visualization with significnace
ggbetweenstats(
data = bottle_join,
x = treatment,
y = final_ta,
type = "parametric", # ANOVA or Kruskal-Wallis
var.equal = TRUE, # ANOVA or Welch ANOVA
plot.type = "box",
pairwise.comparisons = TRUE,
pairwise.display = "significant",
centrality.plotting = FALSE,
bf.message = FALSE
)
Betwen tanks:
Code
# Test tank as a fixed effect WITHIN each treatment
# This directly asks: do my replicate tanks differ?
control_test <- lm(final_ta ~ tank,
data = bottle_join %>% filter(treatment == "ambient"))
summary(control_test)
medium_test <- lm(final_ta ~ tank,
data = bottle_join %>% filter(treatment == "medium"))
summary(medium_test)
high_test <- lm(final_ta ~ tank,
data = bottle_join %>% filter(treatment == "high"))
summary(high_test)No significant differences within replicate treatment tanks.
Treatment stability.
# Create days since experiment start
bottle_join <- bottle_join |>
mutate(day = as.numeric(date - min(date)))
# Look at
treat_time_lmer <- lmer(final_ta ~ treatment * day + (1|tank), data = bottle_join)
summary(treat_time_lmer)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: final_ta ~ treatment * day + (1 | tank)
Data: bottle_join
REML criterion at convergence: 1285.2
Scaled residuals:
Min 1Q Median 3Q Max
-2.6965 -0.4353 0.0021 0.3648 6.1022
Random effects:
Groups Name Variance Std.Dev.
tank (Intercept) 377.2 19.42
Residual 3189.5 56.48
Number of obs: 120, groups: tank, 6
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 2403.1582 23.2878 11.7289 103.194 < 2e-16 ***
treatmentmedium 315.2096 32.5817 11.2657 9.674 8.51e-07 ***
treatmenthigh 705.9799 32.9759 11.7844 21.409 8.49e-11 ***
day 0.8906 0.6008 111.0518 1.482 0.1411
treatmentmedium:day 1.8352 0.8448 111.0322 2.172 0.0319 *
treatmenthigh:day 1.6342 0.8544 111.0096 1.913 0.0584 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) trtmntm trtmnth day trtmntm:
treatmntmdm -0.715
treatmnthgh -0.706 0.505
day -0.713 0.510 0.504
trtmntmdm:d 0.507 -0.704 -0.358 -0.711
trtmnthgh:d 0.502 -0.358 -0.711 -0.703 0.500 So the medium treatment did vary with time a bit.
seacarb
Using pH and TA.
Calculated carbonate chemistry summaries by treatment
Code
# Calculate treatment summary
# Generate carb chem parameters
stats_pHTA <- bottle_join |>
filter(!is.na(glass_p_h), !is.na(final_ta), !is.na(temp)) |>
mutate(
carb =seacarb::carb(
flag = 8,
var1 = glass_p_h,
var2 = final_ta/1e6,
S = coalesce(salinity, 35),
T = temp,
P = 0
),
OmegaAr = carb$OmegaAr,
pCO2 = carb$pCO2,
dic = carb$DIC*1e6,
) |>
filter(!is.na(treatment)) |>
group_by(treatment) |>
summarise(
n = sum(!is.na(final_ta)),
TA_mean = mean(final_ta, na.rm = TRUE),
TA_sd = sd(final_ta, na.rm = TRUE),
TA_sem = TA_sd / sqrt(n),
pH_mean = mean(glass_p_h, na.rm = TRUE),
pH_sd = sd(glass_p_h, na.rm = TRUE),
pH_sem = pH_sd/sqrt(n),
omega_mean = mean(OmegaAr),
omega_sd = sd(OmegaAr),
omega_sem = omega_sd/sqrt(n),
pCO2_mean = mean(pCO2),
pCO2_sd = sd(pCO2),
pCO2_sem = pCO2_sd/sqrt(n),
dic_mean = mean(dic, na.rm = TRUE),
dic_sd = sd(dic, na.rm = TRUE),
dic_sem = dic_sd/sqrt(n),
.groups = "drop"
) |>
arrange(treatment) |>
mutate(TA_delta_tr = TA_mean - lag(TA_mean),
pH_delta_tr = pH_mean - lag(pH_mean),
omega_delta_tr = omega_mean - lag(omega_mean),
pCO2_delta_tr = pCO2_mean - lag(pCO2_mean)
)
print(stats_pHTA, n = Inf, width = Inf)# A tibble: 3 × 21
treatment n TA_mean TA_sd TA_sem pH_mean pH_sd pH_sem omega_mean
<fct> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 ambient 28 2429. 35.4 6.68 8.02 0.0274 0.00517 3.74
2 medium 26 2825. 47.1 9.23 8.07 0.0313 0.00613 4.77
3 high 26 3200. 83.6 16.4 8.11 0.0427 0.00838 5.84
omega_sd omega_sem pCO2_mean pCO2_sd pCO2_sem dic_mean dic_sd dic_sem
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 0.210 0.0398 452. 34.2 6.46 2107. 30.9 5.84
2 0.339 0.0666 457. 35.4 6.94 2431. 30.6 6.01
3 0.480 0.0941 462. 53.8 10.6 2734. 84.0 16.5
TA_delta_tr pH_delta_tr omega_delta_tr pCO2_delta_tr
<dbl> <dbl> <dbl> <dbl>
1 NA NA NA NA
2 396. 0.0525 1.04 4.67
3 375. 0.0445 1.06 4.52Using DIC and TA
Code
# Import DIC data
stats_dicTA <- dic_df_sc |>
mutate(
carb =seacarb::carb(
flag = 15,
var1 = final_ta/1e6,
var2 = dic/1e6,
S = coalesce(salinity, 35),
T = temp,
P = 0
),
OmegaAr = carb$OmegaAr,
pCO2 = carb$pCO2,
pH = carb$pH
) |>
#filter(!is.na(treatment)) |>
group_by(treatment) |>
summarise(
n = sum(!is.na(final_ta)),
TA_mean = mean(final_ta, na.rm = TRUE),
TA_sd = sd(final_ta, na.rm = TRUE),
TA_sem = TA_sd / sqrt(n),
pH_mean = mean(pH, na.rm = TRUE),
pH_sd = sd(pH, na.rm = TRUE),
pH_sem = pH_sd/sqrt(n),
omega_mean = mean(OmegaAr),
omega_sd = sd(OmegaAr),
omega_sem = omega_sd/sqrt(n),
pCO2_mean = mean(pCO2),
pCO2_sd = sd(pCO2),
pCO2_sem = pCO2_sd/sqrt(n),
dic_mean = mean(dic, na.rm = TRUE),
dic_sd = sd(dic, na.rm = TRUE),
dic_sem = dic_sd/sqrt(n),
.groups = "drop"
) |>
arrange(treatment) |>
mutate(TA_delta_tr = TA_mean - lag(TA_mean),
pH_delta_tr = pH_mean - lag(pH_mean),
omega_delta_tr = omega_mean - lag(omega_mean),
pCO2_delta_tr = pCO2_mean - lag(pCO2_mean)
)
print(stats_dicTA, n = Inf, width = Inf)# A tibble: 3 × 21
treatment n TA_mean TA_sd TA_sem pH_mean pH_sd pH_sem omega_mean omega_sd
<fct> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 ambient 10 2422. 46.9 14.8 8.12 0.0942 0.0298 4.49 0.742
2 medium 9 2837. 56.0 18.7 8.15 0.0492 0.0164 5.49 0.537
3 high 9 3326. 415. 138. 8.22 0.0963 0.0321 7.21 1.50
omega_sem pCO2_mean pCO2_sd pCO2_sem dic_mean dic_sd dic_sem TA_delta_tr
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 0.235 341. 88.3 27.9 2029. 75.8 24.0 NA
2 0.179 367. 50.7 16.9 2380. 37.0 12.3 415.
3 0.500 364. 96.6 32.2 2750. 364. 121. 489.
pH_delta_tr omega_delta_tr pCO2_delta_tr
<dbl> <dbl> <dbl>
1 NA NA NA
2 0.0247 0.996 25.9
3 0.0676 1.72 -3.38Compare pH and TA paired samples versus dic and TA pairs
Code
# Add source label and bind
combined_stats <- bind_rows(
stats_pHTA |> mutate(source = "pH_TA"),
stats_dicTA |> mutate(source = "dic_TA")
)
# Pivot to long format
combined_long <- combined_stats |>
dplyr::select(treatment, source, TA_mean, TA_sem, pH_mean, pH_sem, omega_mean, omega_sem, pCO2_mean, pCO2_sem, dic_mean, dic_sem) |>
pivot_longer(
cols = -c(treatment, source),
names_to = c("variable", ".value"),
names_pattern = "(.+)_(mean|sem)"
)
ggplot(combined_long, aes(x = treatment, y = mean, color = source, group = source)) +
geom_point(position = position_dodge(0.3), size = 3) +
geom_errorbar(aes(ymin = mean - sem, ymax = mean + sem),
position = position_dodge(0.3), width = 0.2) +
facet_wrap(~variable, scales = "free_y") +
labs(x = "Treatment", y = NULL, color = "Method") +
theme_classic()
plots
Treatment stability over time:
Code
# Plot TA over time by treatment
ggplot(data = bottle_join |> filter(!is.na(treatment), !is.na(final_ta)),
aes(x = date, y = final_ta, color = treatment, group = treatment)) +
geom_point() +
geom_smooth(data = ~ filter(.x, treatment == "medium"),
method = "lm", se = TRUE) +
theme_bw() +
labs(x = "Date", y = "Total alkalinity (umol/kg)", color = "Treatment") `geom_smooth()` using formula = 'y ~ x'
By tank:
Code
# Plot TA over time by treatment
ggplot(data = bottle_join |> filter(!is.na(treatment), !is.na(final_ta)),
aes(x = date, y = final_ta, color = tank)) +
geom_point() +
geom_smooth(method = "lm", se = TRUE) +
theme_bw() +
labs(x = "Date", y = "Total alkalinity (umol/kg)", color = "Tank") `geom_smooth()` using formula = 'y ~ x'
TA and pH by tank/treatment:
Code
# Boxplots of TA and pH by tank/treatment
ggplot(bottle_long |>
filter(parameter %in% c("final_ta", "glass_p_h"),
!is.na(value)),
aes(x = tank, y = value, fill = treatment)) +
geom_jitter(alpha = 0.3) +
geom_boxplot(alpha = 0.7) +
facet_wrap(~parameter, scales = "free_y") +
theme_classic() +
labs(y = "TA (umol/kg)", fill = "OmegaArg") +
scale_fill_discrete(labels = c("3.73", "4.79", "5.96"))
Growth
Change in planar area of the recruits was tracked with ImageJ. A 0.2 mm ruler tick was used to set the scale. The same 1 mm tick mark was measured three times in each photo to collect method precision estimates. Each coral was traced by hand three times and the area and perimeter recorded. Size at each timepoint is the mean of the three replicate area measurements.
ImageJ method precision and accuracy
Standard was measured with straight line tool. Convert linear uncertainty to area uncertainty
Code
relError <- grData |>
filter(!is.na(standard_mm)) |>
summarise(stand_mean = mean(standard_mm),
stand_sd = sd(standard_mm),
n = n()) |>
mutate(linSD= stand_sd/stand_mean) |> # Calculate linear uncertainty
mutate(areaSD = linSD * 2) # Transform to area uncertaintyThe standard mean = 1.01, SD = 0.01
The fractional uncertainty in the length measurement:
Relative linear uncertainty = SD of standard / mean length of standard * 100 = 1.04%
The relative uncertainty in area is twice the relative uncertainty in length, so: area uncertainty = 2 x 1.04 = 2.07%
Plots
Planar area
Growth rate
Growth rate (change in relative area per day)
Code
grSumm <- grSumm |>
group_by(coral) |>
arrange(date) |>
mutate(
initial_area = first(area_mean),
days = as.numeric(date - first(date)),
growth_rate = (area_mean - initial_area) / initial_area / days
)
grSumm |>
filter(!is.nan(growth_rate)) |>
group_by(species, treatment, date) |>
summarise(
mean_rate = mean(growth_rate, na.rm = TRUE),
se_rate = sd(growth_rate, na.rm = TRUE) / sqrt(sum(!is.na(growth_rate))),
.groups = "drop"
) |>
ggplot(aes(x = date, y = mean_rate, color = treatment, group = treatment)) +
geom_line() +
geom_point() +
geom_errorbar(aes(ymin = mean_rate - se_rate, ymax = mean_rate + se_rate), width = 0.2) +
facet_wrap(~species) +
labs(y = "Growth rate (relative area / day)", x = "Date") +
theme_bw()
Percent increase
Percent change in area across timepoints
Code
# # Percent increase across timepoints
grSumm_long <- grSumm |>
dplyr::select(coral, tp, area_mean, perim_mean, species, days, tank, initial_area) |>
dplyr::arrange(coral, tp) |>
dplyr::group_by(coral) |>
# baseline
dplyr::mutate(area_tp0 = area_mean[tp == 0]) |>
# regular percent increase from baseline
dplyr::mutate(
percent_increase_raw = (area_mean - area_tp0) / area_tp0 * 100
) |>
# force it to never decrease (carry forward previous max)
dplyr::mutate(
percent_increase = cummax(percent_increase_raw)
) |>
dplyr::ungroup()
grSumm_long <- left_join(
grSumm_long,
coral_metadata |>
dplyr::select(coral, treatment),
by = "coral"
)
grSumm_long <- grSumm_long |>
mutate(treatment = as.character(treatment), # make sure we can assign new value
treatment = if_else(is.na(treatment), "ambient", treatment), # replace NAs
treatment = factor(treatment, levels = c("ambient", "medium", "high")))Code
grSumm_long |>
group_by(treatment, tp, species) |>
summarise(
mean_perc = mean(percent_increase, na.rm = TRUE),
se_perc = sd(percent_increase, na.rm = TRUE)/sqrt(n()),
.groups = "drop" ) |>
ggplot(aes(x = tp, y = mean_perc, color = treatment)) +
geom_line(aes(group = treatment), linetype = 1) +
geom_point(size = 3) +
geom_errorbar(aes(ymin = mean_perc - se_perc, ymax = mean_perc + se_perc), width = 0.2) +
theme_bw() +
facet_wrap(~species) +
labs(y = "Percent increase", x = "Timepoint")
Code
# by tank
grSumm_long |>
dplyr::filter(species == "pstr",
!is.na(tank),
!is.na(species),
!is.na(treatment),
!is.na(percent_increase),) |>
group_by(treatment, tp, tank) |>
summarise(
mean_perc = mean(percent_increase, na.rm = TRUE),
se_perc = sd(percent_increase, na.rm = TRUE)/sqrt(n()),
.groups = "drop" ) |>
ggplot(aes(x = tp, y = mean_perc, color = treatment)) +
geom_line(aes(group = treatment), linetype = 1) +
geom_point(size = 3) +
geom_errorbar(aes(ymin = mean_perc - se_perc, ymax = mean_perc + se_perc), width = 0.2) +
theme_bw() +
facet_wrap(~tank) +
labs(y = "Percent increase", x = "Timepoint", title = "PSTR percent increase in area by tank")
Code
grSumm_long |>
dplyr::filter(species == "dlab",
!is.na(tank),
!is.na(species),
!is.na(treatment),
!is.na(percent_increase),) |>
group_by(treatment, tp, tank) |>
summarise(
mean_perc = mean(percent_increase, na.rm = TRUE),
se_perc = sd(percent_increase, na.rm = TRUE)/sqrt(n()),
.groups = "drop" ) |>
ggplot(aes(x = tp, y = mean_perc, color = treatment)) +
geom_line(aes(group = treatment), linetype = 1) +
geom_point(size = 3) +
geom_errorbar(aes(ymin = mean_perc - se_perc, ymax = mean_perc + se_perc), width = 0.2) +
theme_bw() +
facet_wrap(~tank) +
labs(y = "Percent increase", x = "Timepoint", title = "DLAB percent increase in area by tank")
Table Summary:
Code
grSumm_long |>
group_by(treatment, species, tp) |>
summarise(mean_perc_increase = mean(percent_increase, na.rm = TRUE),
sd_perc_increase = sd(percent_increase, na.rm = TRUE),
n = n(),
.groups = "drop_last" # Removes message from output
)# A tibble: 24 × 6
# Groups: treatment, species [6]
treatment species tp mean_perc_increase sd_perc_increase n
<fct> <chr> <int> <dbl> <dbl> <int>
1 ambient dlab 0 0 0 16
2 ambient dlab 4 4.57 4.93 16
3 ambient dlab 5 7.76 5.87 16
4 ambient dlab 7 14.0 13.3 16
5 ambient pstr 0 0 0 15
6 ambient pstr 4 24.3 23.3 15
7 ambient pstr 5 30.2 21.7 15
8 ambient pstr 7 50.3 31.5 15
9 medium dlab 0 0 0 14
10 medium dlab 4 8.19 9.92 14
# ℹ 14 more rowsraw area by days
Code
grSumm |>
group_by(species, treatment, date) |>
summarise(
mean_rate = mean(area_mean, na.rm = TRUE),
se_rate = sd(area_mean, na.rm = TRUE) / sqrt(sum(!is.na(area_mean))),
.groups = "drop"
) |>
ggplot(aes(x = date, y = mean_rate, color = treatment, group = treatment)) +
geom_line() +
geom_point() +
geom_errorbar(aes(ymin = mean_rate - se_rate, ymax = mean_rate + se_rate), width = 0.2) +
facet_wrap(~species) +
labs(y = "Area (mm^2)", x = "Date") +
theme_bw()
Stats
PSTR.
LMM
Run separate linear models with growth rate and raw area as response variables. Treatment, tank, days, and their interaction as fixed effect, and tank and coral as nested random effects.
Code
# using area_mean
area_lmer <- lmer(area_mean ~ treatment * days + initial_area + (1 | tank/coral),
data = grSumm %>% filter(species == "pstr"))
# using growth rate
rate_lmer <- lmer(growth_rate ~ treatment * days + (1 | tank/coral),
data = grSumm %>% filter(species == "pstr"))
# using percent increase
perc_lmer <- lmer(percent_increase ~ treatment * days + (1 | tank/coral),
data = grSumm_long %>% filter(species == "pstr"))growth rate by treatment
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: growth_rate ~ treatment * days + (1 | tank/coral)
Data: grSumm %>% filter(species == "pstr")
REML criterion at convergence: -789.8
Scaled residuals:
Min 1Q Median 3Q Max
-2.58744 -0.40687 -0.01958 0.35824 3.11222
Random effects:
Groups Name Variance Std.Dev.
coral:tank (Intercept) 3.984e-05 0.006312
tank (Intercept) 0.000e+00 0.000000
Residual 1.341e-05 0.003661
Number of obs: 114, groups: coral:tank, 38; tank, 6
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 5.993e-03 2.783e-03 1.064e+02 2.154 0.0335 *
treatmentmedium -1.429e-03 4.400e-03 1.064e+02 -0.325 0.7459
treatmenthigh 1.857e-03 4.084e-03 1.064e+02 0.455 0.6502
days 6.172e-05 5.811e-05 7.300e+01 1.062 0.2917
treatmentmedium:days 1.448e-04 9.188e-05 7.300e+01 1.577 0.1192
treatmenthigh:days 1.606e-05 8.528e-05 7.300e+01 0.188 0.8512
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) trtmntm trtmnth days trtmntm:
treatmntmdm -0.632
treatmnthgh -0.681 0.431
days -0.787 0.497 0.536
trtmntmdm:d 0.497 -0.787 -0.339 -0.632
trtmnthgh:d 0.536 -0.339 -0.787 -0.681 0.431
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')percent increase in area by treatment
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: percent_increase ~ treatment * days + (1 | tank/coral)
Data: grSumm_long %>% filter(species == "pstr")
REML criterion at convergence: 1326.3
Scaled residuals:
Min 1Q Median 3Q Max
-2.8966 -0.5544 -0.0592 0.4778 3.3108
Random effects:
Groups Name Variance Std.Dev.
coral:tank (Intercept) 2.652e+02 16.284236
tank (Intercept) 5.477e-08 0.000234
Residual 2.271e+02 15.070997
Number of obs: 154, groups: coral:tank, 40; tank, 6
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) -1.09362 5.52799 70.11343 -0.198 0.8437
treatmentmedium -2.52891 8.31871 70.80220 -0.304 0.7620
treatmenthigh -0.05536 8.11287 70.11344 -0.007 0.9946
days 0.96599 0.10675 112.10155 9.049 5.12e-15 ***
treatmentmedium:days 0.47286 0.16641 116.31876 2.842 0.0053 **
treatmenthigh:days 0.19187 0.15667 112.10155 1.225 0.2233
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) trtmntm trtmnth days trtmntm:
treatmntmdm -0.665
treatmnthgh -0.681 0.453
days -0.546 0.363 0.372
trtmntmdm:d 0.350 -0.527 -0.238 -0.642
trtmnthgh:d 0.372 -0.247 -0.546 -0.681 0.437
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')raw area increase by treatment
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: area_mean ~ treatment * days + initial_area + (1 | tank/coral)
Data: grSumm %>% filter(species == "pstr")
REML criterion at convergence: -5.9
Scaled residuals:
Min 1Q Median 3Q Max
-3.1278 -0.4478 -0.0631 0.3570 5.6006
Random effects:
Groups Name Variance Std.Dev.
coral:tank (Intercept) 0.01761 0.1327
tank (Intercept) 0.00000 0.0000
Residual 0.03160 0.1778
Number of obs: 154, groups: coral:tank, 40; tank, 6
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) -2.553e-01 7.672e-02 5.663e+01 -3.328 0.00154 **
treatmentmedium -9.466e-02 8.361e-02 8.865e+01 -1.132 0.26063
treatmenthigh 2.198e-02 8.009e-02 9.019e+01 0.274 0.78434
days 6.855e-03 1.259e-03 1.117e+02 5.445 3.12e-07 ***
initial_area 1.319e+00 7.148e-02 3.636e+01 18.455 < 2e-16 ***
treatmentmedium:days 8.617e-03 1.953e-03 1.164e+02 4.412 2.30e-05 ***
treatmenthigh:days 2.354e-04 1.848e-03 1.117e+02 0.127 0.89887
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) trtmntm trtmnth days intl_r trtmntm:
treatmntmdm -0.327
treatmnthgh -0.529 0.430
days -0.464 0.425 0.444
initial_are -0.704 -0.192 0.065 0.000
trtmntmdm:d 0.329 -0.614 -0.289 -0.645 -0.043
trtmnthgh:d 0.316 -0.290 -0.652 -0.681 0.000 0.439
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')Percent area increase and raw area increase models both suggest medium alkalinity corals grew significantly faster than controls (raw area produces higher significance). High alkalinity corals did not differ significantly from controls, suggesting a non-linear response to alkalinity-enhancement.
GAM
Run GAM to piece apart further – are growth trajectories non-linear across treatments or time? Move forward with raw area as response
Model: Generalized Additive Model with separate smooth growth curve per treatment, coral and tank as random effects
Family: gaussian
Link function: identity
Formula:
area_mean ~ treatment + initial_area + s(days, by = treatment,
k = 4) + s(tank, bs = "re") + s(coral, bs = "re")
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.06669 0.06770 -0.985 0.327
treatmentmedium 0.14184 0.06576 2.157 0.033 *
treatmenthigh 0.02872 0.06050 0.475 0.636
initial_area 1.32265 0.07119 18.580 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(days):treatmentambient 1.323e+00 1.558 21.591 2.13e-06 ***
s(days):treatmentmedium 2.067e+00 2.390 48.680 < 2e-16 ***
s(days):treatmenthigh 1.000e+00 1.000 28.800 1.05e-06 ***
s(tank) 1.724e-05 3.000 0.000 0.528
s(coral) 2.463e+01 36.000 2.276 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
R-sq.(adj) = 0.918 Deviance explained = 93.5%
-REML = -13.222 Scale est. = 0.030156 n = 154
Results. Model fit explains 94% of deviance. All three treatments show significant growth over time. Control and high alkalinity grow linearly. Medium alkalinity shows a slightly accelerating, exponential-like trajectory. Tank variation is negligible — tanks within treatments were consistent. Individual coral identity drives most of the variation in growth.
Medium alkalinity enhanced growth most clearly, with an accelerating trajectory over time. High alkalinity did not differ meaningfully from control.
DLAB.
LMM
Model: Linear mixed effects model with treatment × days interaction as fixed effects, coral nested in tank as random effects
raw area by treatment
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: area_mean ~ treatment * days + (1 | tank/coral)
Data: grSumm %>% filter(species == "dlab")
REML criterion at convergence: 613.9
Scaled residuals:
Min 1Q Median 3Q Max
-3.7365 -0.3203 -0.0630 0.3282 4.1581
Random effects:
Groups Name Variance Std.Dev.
coral:tank (Intercept) 44.1351 6.6434
tank (Intercept) 22.0920 4.7002
Residual 0.3529 0.5941
Number of obs: 177, groups: coral:tank, 45; tank, 6
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 8.910171 3.737478 2.949838 2.384 0.09874 .
treatmentmedium -1.283430 5.314109 3.019236 -0.242 0.82463
treatmenthigh -0.537525 5.290318 2.968029 -0.102 0.92554
days 0.009500 0.004131 129.003794 2.300 0.02306 *
treatmentmedium:days 0.013070 0.006106 129.005792 2.140 0.03420 *
treatmenthigh:days 0.019832 0.005897 129.001594 3.363 0.00101 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) trtmntm trtmnth days trtmntm:
treatmntmdm -0.703
treatmnthgh -0.706 0.497
days -0.031 0.022 0.022
trtmntmdm:d 0.021 -0.032 -0.015 -0.676
trtmnthgh:d 0.022 -0.015 -0.031 -0.701 0.474 growth rate by treatment
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: growth_rate ~ treatment * days + (1 | tank/coral)
Data: grSumm %>% filter(species == "dlab")
REML criterion at convergence: -1102
Scaled residuals:
Min 1Q Median 3Q Max
-2.65035 -0.41249 -0.04961 0.36078 3.08746
Random effects:
Groups Name Variance Std.Dev.
coral:tank (Intercept) 1.332e-05 0.003650
tank (Intercept) 1.077e-06 0.001038
Residual 2.822e-06 0.001680
Number of obs: 132, groups: coral:tank, 45; tank, 6
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 1.450e-04 1.557e-03 6.916e+00 0.093 0.928
treatmentmedium 1.674e-03 2.268e-03 7.918e+00 0.738 0.482
treatmenthigh 5.691e-03 2.219e-03 7.301e+00 2.564 0.036 *
days 4.332e-05 2.640e-05 8.393e+01 1.641 0.105
treatmentmedium:days -2.378e-05 3.930e-05 8.407e+01 -0.605 0.547
treatmenthigh:days -7.183e-05 3.752e-05 8.382e+01 -1.914 0.059 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) trtmntm trtmnth days trtmntm:
treatmntmdm -0.687
treatmnthgh -0.702 0.482
days -0.633 0.435 0.444
trtmntmdm:d 0.425 -0.642 -0.299 -0.672
trtmnthgh:d 0.446 -0.306 -0.634 -0.704 0.473 Results
All treatments grow significantly over time At day 0, no significant size differences between treatments — good starting conditions Control grows at 0.0095 standardized units/day Medium alkalinity grows significantly faster than control (p = 0.034) High alkalinity grows significantly faster than control (p = 0.001) — stronger effect than medium Combined growth rates: control 0.0095, medium 0.0226, high 0.0293/day
Unlike pstr, both alkalinity treatments significantly enhanced dlab growth, with high alkalinity showing the strongest effect — suggesting dlab may respond more positively to elevated alkalinity across the full range tested.
GAM
Model: Generalized Additive Model with separate smooth growth curve per treatment, coral and tank as random effects
Family: gaussian
Link function: identity
Formula:
area_mean ~ treatment + s(days, by = treatment, k = 4) + s(tank,
bs = "re") + s(coral, bs = "re")
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 9.1760 3.7384 2.455 0.0154 *
treatmentmedium -0.9124 5.3152 -0.172 0.8640
treatmenthigh 0.0149 5.2915 0.003 0.9978
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(days):treatmentambient 1.007 1.014 5.42 0.02149 *
s(days):treatmentmedium 2.065 2.394 12.24 4.94e-06 ***
s(days):treatmenthigh 1.000 1.000 50.17 < 2e-16 ***
s(tank) 2.354 3.000 724960.04 0.00362 **
s(coral) 39.568 42.000 9730.10 0.24531
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
R-sq.(adj) = 0.994 Deviance explained = 99.5%
-REML = 297.16 Scale est. = 0.34171 n = 177
Results.
Model fexplains 99.5% of deviance All three treatments show significant growth over time Control and high alkalinity grow linearly (edf ≈ 1) Medium alkalinity shows a slightly non-linear, accelerating trajectory (edf = 2.07) — similar pattern to pstr Tank variation is significant (p = 0.004) — some between-tank differences worth noting Coral identity not significant here (p = 0.245) — unlike pstr, individual coral variation is less important
Both alkalinity treatments enhanced dlab growth, with high alkalinity showing the strongest and most linear response. Medium again shows an accelerating curve over time, consistent with what was seen in pstr.
Photophysiology
DLAB
Stats
ANOVA and Tukey
Code
dlabYIIaov <- aov(y ~ treatment,
data = pam_df |>
filter(species == "dlab")
)
summary(dlabYIIaov) Df Sum Sq Mean Sq F value Pr(>F)
treatment 2 0.06231 0.031157 23.33 7.8e-10 ***
Residuals 200 0.26706 0.001335
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
21 observations deleted due to missingnessCode
dlabYIItukey <- TukeyHSD(dlabYIIaov)
print(dlabYIItukey) Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = y ~ treatment, data = filter(pam_df, species == "dlab"))
$treatment
diff lwr upr p adj
medium-ambient -0.02289240 -0.03779446 -0.007990342 0.0010548
high-ambient -0.04196878 -0.05651098 -0.027426569 0.0000000
high-medium -0.01907638 -0.03422796 -0.003924793 0.0092468Linear model
Code
dlabYIIlm <- lm(y ~ treatment + tank,
data = pam_df,
subset = species == "dlab")
anova(dlabYIIlm)Analysis of Variance Table
Response: y
Df Sum Sq Mean Sq F value Pr(>F)
treatment 2 0.062313 0.0311567 24.9576 2.182e-10 ***
tank 3 0.021124 0.0070415 5.6405 0.000996 ***
Residuals 197 0.245932 0.0012484
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Code
emmeans(dlabYIIlm, pairwise ~ tank | treatment, adjust = "tukey")NOTE: A nesting structure was detected in the fitted model:
tank %in% treatment$emmeans
treatment = ambient:
tank emmean SE df lower.CL upper.CL
B1 0.533 0.00721 197 0.519 0.547
D1 0.548 0.00505 197 0.539 0.558
treatment = medium:
tank emmean SE df lower.CL upper.CL
B2 0.521 0.00625 197 0.508 0.533
D2 0.520 0.00645 197 0.508 0.533
treatment = high:
tank emmean SE df lower.CL upper.CL
B3 0.521 0.00668 197 0.507 0.534
D3 0.488 0.00559 197 0.477 0.499
Confidence level used: 0.95
$contrasts
treatment = ambient:
contrast estimate SE df t.ratio p.value
B1 - D1 -0.01534 0.00880 197 -1.743 0.0829
treatment = high:
contrast estimate SE df t.ratio p.value
B3 - D3 0.03244 0.00871 197 3.726 0.0003
treatment = medium:
contrast estimate SE df t.ratio p.value
B2 - D2 0.00045 0.00898 197 0.050 0.9601PSTR
Stats
ANOVA and Tukey
Code
pstrYIIaov <- aov(y ~ treatment,
data = pam_df |>
filter(species == "pstr")
)
summary(pstrYIIaov) Df Sum Sq Mean Sq F value Pr(>F)
treatment 2 0.00385 0.001927 1.272 0.283
Residuals 189 0.28629 0.001515 Code
pstrYIItukey <- TukeyHSD(pstrYIIaov)
print(pstrYIItukey) Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = y ~ treatment, data = filter(pam_df, species == "pstr"))
$treatment
diff lwr upr p adj
medium-ambient -0.011083590 -0.02767497 0.005507787 0.2576123
high-ambient -0.005906667 -0.02148746 0.009674127 0.6437768
high-medium 0.005176923 -0.01192904 0.022282883 0.7549712Linear model
Code
pstrYIIlm <- lm(y ~ treatment + tank,
data = pam_df,
subset = species == "pstr")
anova(pstrYIIlm)Analysis of Variance Table
Response: y
Df Sum Sq Mean Sq F value Pr(>F)
treatment 2 0.003853 0.0019265 1.4636 0.2341
tank 3 0.041468 0.0138225 10.5014 2.052e-06 ***
Residuals 186 0.244823 0.0013163
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Code
emmeans(pstrYIIlm, pairwise ~ tank | treatment, adjust = "tukey")NOTE: A nesting structure was detected in the fitted model:
tank %in% treatment$emmeans
treatment = ambient:
tank emmean SE df lower.CL upper.CL
B1 0.560 0.00541 186 0.549 0.571
D1 0.524 0.00662 186 0.511 0.537
treatment = medium:
tank emmean SE df lower.CL upper.CL
B2 0.542 0.00698 186 0.528 0.556
D2 0.527 0.00726 186 0.512 0.541
treatment = high:
tank emmean SE df lower.CL upper.CL
B3 0.556 0.00662 186 0.543 0.569
D3 0.526 0.00613 186 0.513 0.538
Confidence level used: 0.95
$contrasts
treatment = ambient:
contrast estimate SE df t.ratio p.value
B1 - D1 0.0361 0.00855 186 4.218 <0.0001
treatment = high:
contrast estimate SE df t.ratio p.value
B3 - D3 0.0305 0.00903 186 3.381 0.0009
treatment = medium:
contrast estimate SE df t.ratio p.value
B2 - D2 0.0152 0.01010 186 1.512 0.1324Health over time
Code
survivor_df<- left_join(
grSumm,
coral_metadata |>
dplyr::select(coral),
by = "coral"
)
survivor_df <- survivor_df %>%
mutate(treatment = as.character(treatment), # make sure we can assign new value
treatment = if_else(is.na(treatment), "ambient", treatment), # replace NAs
treatment = factor(treatment, levels = c("ambient", "medium", "high")))
survivor_df <- survivor_df |>
group_by(tp, treatment, health) |>
summarise(n = n(), .groups = "drop") |>
group_by(tp, treatment) |>
mutate(prop = n / sum(n))
survivor_df <- survivor_df |>
mutate(health = na_if(health, ""))
ggplot(survivor_df |> filter(!is.na(health)),
aes(x = tp, y = prop, fill = health)) +
geom_bar(stat = "identity") +
scale_fill_discrete(labels = c("dead", "healthy", "tissue recession")) +
facet_wrap(~treatment) +
theme_classic()