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.

# 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

# 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)
tidy(treat_lmer, effects = "fixed") |> knitr::kable(digits = 3)

effect	term	estimate	std.error	statistic	df	p.value
fixed	(Intercept)	2427.782	34.252	70.881	2.945	0.000
fixed	treatmentmedium	363.295	48.541	7.484	2.970	0.005
fixed	treatmenthigh	781.937	48.548	16.106	2.971	0.001

report(treat_lmer)

We fitted a linear mixed model (estimated using REML and nloptwrap optimizer)
to predict final_ta with treatment (formula: final_ta ~ treatment). The model
included tank as random effect (formula: ~1 | tank). The model's total
explanatory power is substantial (conditional R2 = 0.86) and the part related
to the fixed effects alone (marginal R2) is of 0.85. The model's intercept,
corresponding to treatment = ambient, is at 2427.78 (95% CI [2359.94, 2495.62],
t(116) = 70.88, p < .001). Within this model:

  - The effect of treatment [medium] is statistically significant and positive
(beta = 363.29, 95% CI [267.15, 459.44], t(116) = 7.48, p < .001; Std. beta =
1.05, 95% CI [0.77, 1.33])
  - The effect of treatment [high] is statistically significant and positive
(beta = 781.94, 95% CI [685.78, 878.09], t(116) = 16.11, p < .001; Std. beta =
2.26, 95% CI [1.98, 2.54])

Standardized parameters were obtained by fitting the model on a standardized
version of the dataset. 95% Confidence Intervals (CIs) and p-values were
computed using a Wald t-distribution approximation.

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:

# Residuals vs fitted 
plot(treat_lmer)

# 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.

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)
tidy(treat_lmer, effects = "fixed") |> knitr::kable(digits = 3)

effect	term	estimate	std.error	statistic	df	p.value
fixed	(Intercept)	2427.787	16.219	149.687	2.879	0.000
fixed	treatmentmedium	363.289	22.993	15.800	2.907	0.001
fixed	treatmenthigh	749.979	23.060	32.523	2.939	0.000

report(treat_lmer)

We fitted a linear mixed model (estimated using REML and nloptwrap optimizer)
to predict final_ta with treatment (formula: final_ta ~ treatment). The model
included tank as random effect (formula: ~1 | tank). The model's total
explanatory power is substantial (conditional R2 = 0.96) and the part related
to the fixed effects alone (marginal R2) is of 0.95. The model's intercept,
corresponding to treatment = ambient, is at 2427.79 (95% CI [2395.66, 2459.91],
t(115) = 149.69, p < .001). Within this model:

  - The effect of treatment [medium] is statistically significant and positive
(beta = 363.29, 95% CI [317.74, 408.83], t(115) = 15.80, p < .001; Std. beta =
1.16, 95% CI [1.01, 1.30])
  - The effect of treatment [high] is statistically significant and positive
(beta = 749.98, 95% CI [704.30, 795.66], t(115) = 32.52, p < .001; Std. beta =
2.39, 95% CI [2.24, 2.54])

Standardized parameters were obtained by fitting the model on a standardized
version of the dataset. 95% Confidence Intervals (CIs) and p-values were
computed using a Wald t-distribution approximation.

# 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

treat_aov <- aov(final_ta ~ treatment, data = bottle_join)
broom::tidy(treat_aov) |> knitr::kable(digits = 3)

term	df	sumsq	meansq	statistic	p.value
treatment	2	11203938.8	5601969.42	1272.698	0
Residuals	117	514993.1	4401.65	NA	NA

# Post-hoc comparison
treat_tukey <- glht(treat_aov,
                    linfct = mcp(treatment = "Tukey"))

broom::tidy(treat_tukey) |> knitr::kable(digits = 3)

term	contrast	null.value	estimate	std.error	statistic	adj.p.value
treatment	medium - ambient	0	363.223	14.744	24.635	0
treatment	high - ambient	0	748.694	14.840	50.452	0
treatment	high - medium	0	385.471	14.930	25.819	0

cld(treat_tukey) |> broom::tidy() |> knitr::kable(digits = 3)

treatment	letters
ambient	a
medium	b
high	c

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.16

# 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:

# 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"))
broom::tidy(control_test) |> knitr::kable(digits = 3)

medium_test <- lm(final_ta ~ tank, 
                  data = bottle_join %>% filter(treatment == "medium"))
broom::tidy(medium_test) |> knitr::kable(digits = 3)

high_test <- lm(final_ta ~ tank, 
                data = bottle_join %>% filter(treatment == "high"))
broom::tidy(high_test) |> knitr::kable(digits = 3)

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)))

# LMER  -- "pulls out linear drift, so if last tp sig higher than first, will see linear effect
treat_time_lmer <- lmer(final_ta ~ treatment * day + (1|tank), data = bottle_join)
tidy(treat_time_lmer, effects = "fixed") |> 
  knitr::kable(digits = 3, caption = "Fixed Effects: TA ~ Treatment x Day")

Fixed Effects: TA ~ Treatment x Day

effect	term	estimate	std.error	statistic	df	p.value
fixed	(Intercept)	2403.158	23.288	103.194	11.729	0.000
fixed	treatmentmedium	315.210	32.582	9.674	11.266	0.000
fixed	treatmenthigh	705.980	32.976	21.409	11.784	0.000
fixed	day	0.891	0.601	1.482	111.052	0.141
fixed	treatmentmedium:day	1.835	0.845	2.172	111.032	0.032
fixed	treatmenthigh:day	1.634	0.854	1.913	111.010	0.058

# GAMM -- might better handle variability around the treatments
gamTreat <- gam(final_ta ~ treatment + s(day, by = treatment, k = 4) +
              s(tank, bs = "re"),
            data = bottle_join,
            method = "REML")

draw(gamTreat)  # plots smooths cleanly

Linear model shows the medium treatment did vary with time a bit. GAMM model shows significant variation over time in all treatments except ambient.

seacarb

Using pH and TA.

Calculated carbonate chemistry summaries by treatment

# 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)
         )


kable(stats_pHTA, digits = 3, caption = "Carb Chem Summary using pH and TA")

Carb Chem Summary using pH and TA

treatment	n	TA_mean	TA_sd	TA_sem	pH_mean	pH_sd	pH_sem	omega_mean	omega_sd	omega_sem	pCO2_mean	pCO2_sd	pCO2_sem	dic_mean	dic_sd	dic_sem	TA_delta_tr	pH_delta_tr	omega_delta_tr	pCO2_delta_tr
ambient	28	2429.272	35.353	6.681	8.017	0.027	0.005	3.738	0.210	0.040	452.437	34.205	6.464	2107.151	30.916	5.843	NA	NA	NA	NA
medium	26	2825.317	47.067	9.231	8.070	0.031	0.006	4.775	0.339	0.067	457.111	35.380	6.939	2430.640	30.620	6.005	396.045	0.053	1.037	4.673
high	26	3200.038	83.636	16.402	8.114	0.043	0.008	5.835	0.480	0.094	461.631	53.811	10.553	2734.214	83.992	16.472	374.721	0.044	1.060	4.520

Using DIC and TA

# 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)
         )

kable(stats_dicTA, digits = 3, caption = "Carb Chem Summary using DIC and TA")

Carb Chem Summary using DIC and TA

treatment	n	TA_mean	TA_sd	TA_sem	pH_mean	pH_sd	pH_sem	omega_mean	omega_sd	omega_sem	pCO2_mean	pCO2_sd	pCO2_sem	dic_mean	dic_sd	dic_sem	TA_delta_tr	pH_delta_tr	omega_delta_tr	pCO2_delta_tr
ambient	10	2422.174	46.862	14.819	8.123	0.094	0.030	4.494	0.742	0.235	341.258	88.340	27.936	2028.838	75.784	23.965	NA	NA	NA	NA
medium	9	2836.976	56.004	18.668	8.148	0.049	0.016	5.490	0.537	0.179	367.159	50.746	16.915	2379.524	37.041	12.347	414.802	0.025	0.996	25.901
high	9	3326.116	415.161	138.387	8.215	0.096	0.032	7.213	1.500	0.500	363.779	96.557	32.186	2749.641	363.548	121.183	489.141	0.068	1.723	-3.380

Compare pH and TA paired samples versus dic and TA pairs

# 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:

# 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:

# 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:

# 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

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 uncertainty

The 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)

# Calculate growth rates per interval
grSumm <- grSumm |>
  group_by(coral) |>
  arrange(date) |>
  mutate(
    initial_area = first(area_mean),
    days = as.numeric(date - first(date)),
    growth_rate = (area_mean - lag(area_mean)) / initial_area / (days - lag(days))
  )


# Plot
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()

Warning: Removed 6 rows containing missing values or values outside the scale range
(`geom_line()`).

Warning: Removed 6 rows containing missing values or values outside the scale range
(`geom_point()`).

# Cumulative summary
rateCumulative <- grSumm |> 
  filter(tp == 7) |> 
  mutate(growth_rate_cumulative = (area_mean - initial_area) / initial_area / days) |> 
  group_by(treatment, species) |> 
  summarise(
    n = n(),
    mean_rate = mean(growth_rate_cumulative, na.rm = TRUE),
    sd_rate = sd(growth_rate_cumulative, na.rm = TRUE),
    se_rate = sd_rate/sqrt(n),
    .groups = "drop"
  )

rateTreatDiff <- rateCumulative |> 
  group_by(species) |> 
  mutate(
    amb_rate = mean_rate[treatment == "ambient"], 
    pct_change = (mean_rate - amb_rate) / amb_rate * 100
  )

# Define treatment effects 
pstr_med_pct <- rateTreatDiff |> filter(species == "pstr", treatment == "medium") |> pull(pct_change) |> round(1)
pstr_high_pct <-  rateTreatDiff |> filter(species == "pstr", treatment == "high") |> pull(pct_change) |> round(1)

dlab_med_pct <- rateTreatDiff |> filter(species == "dlab", treatment == "medium") |> pull(pct_change) |> round(1)
dlab_high_pct <- rateTreatDiff |> filter(species == "dlab", treatment == "high") |> pull(pct_change) |> round(1)

Mean growth rates per treatment (mm^2/day)

Mean growth rates per treatment

treatment	species	n	mean_rate	sd_rate	se_rate
ambient	dlab	16	0.0026	0.0028	0.0007
ambient	pstr	15	0.0102	0.0065	0.0017
medium	dlab	14	0.0041	0.0052	0.0014
medium	pstr	12	0.0152	0.0056	0.0016
high	dlab	15	0.0045	0.0038	0.0010
high	pstr	13	0.0118	0.0071	0.0020

Treatment effect summary on growth rates:
PSTR growth rate increased by 49% in medium treatment and 15.3% in high.
DLAB growth rate increawsed by 58.1% in medium treatment and 73% in high.

Percent increase

Percent change in area across timepoints

# # 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")))

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")

# 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")

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")

percCumulative <- grSumm_long |> 
  filter(tp == 7) |> 
  group_by(treatment, species) |> 
  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
  )

percTreatDiff <- percCumulative |> 
  group_by(species) |> 
  mutate(
    amb_perc = mean_perc_increase[treatment == "ambient"], 
    point_diff = mean_perc_increase - amb_perc
  )

# Define overall percent increase in growth by treatment
pstr_med_perc_area <- percCumulative |> filter(species == "pstr", treatment == "medium") |> pull(mean_perc_increase) |> round(1)
pstr_high_perc_area <- percCumulative |> filter(species == "pstr", treatment == "high") |> pull(mean_perc_increase) |> round(1)

dlab_med_perc_area <- percCumulative |> filter(species == "dlab", treatment == "medium") |> pull(mean_perc_increase) |> round(1)
dlab_high_perc_area <- percCumulative |> filter(species == "dlab", treatment == "high") |> pull(mean_perc_increase) |> round(1)

# Define treatment differences in pp, percentage points
pstr_med_pp <- percTreatDiff |> filter(species == "pstr", treatment == "medium") |> pull(point_diff) |> round(1)
pstr_high_pp <-  percTreatDiff |> filter(species == "pstr", treatment == "high") |> pull(point_diff) |> round(1)

dlab_med_pp <- percTreatDiff |> filter(species == "dlab", treatment == "medium") |> pull(point_diff) |> round(1)
dlab_high_pp <- percTreatDiff |> filter(species == "dlab", treatment == "high") |> pull(point_diff) |> round(1)

Overall, corals exposed to AE treatments showed increased area gain relative to ambient controls across both species. In P. strigosa, medium and high AE treatments resulted in mean percent area increases of 74.7% and 57.9%, representing 24.4 and 7.6 percentage points above ambient controls, respectively. D. labyrinthiformis similarly showed increased area gain under AE treatment, with medium and high treatments reaching 20.9% and 23% area increase, or 6.9 and 9 percentage points above controls, respectively.

Seems like PSTR were growing faster than the DLABs in general, but the treatment effect was a bit stronger and more linear in DLAB.

raw area by days

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.

# 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

Fixed Effects: Growth rate ~ Treatment x Days

effect	term	estimate	std.error	statistic	df	p.value
fixed	(Intercept)	-0.014	0.009	-1.676	108	0.097
fixed	treatmentmedium	0.002	0.014	0.158	108	0.874
fixed	treatmenthigh	0.017	0.013	1.350	108	0.180
fixed	days	0.001	0.000	3.081	108	0.003
fixed	treatmentmedium:days	0.000	0.000	0.315	108	0.754
fixed	treatmenthigh:days	0.000	0.000	-1.277	108	0.204

Corals were all growing over time, but treatment had no effect when looking at rates alone.

percent increase in area by treatment

Fixed Effects: Percent increase ~ Treatment x Days

effect	term	estimate	std.error	statistic	df	p.value
fixed	(Intercept)	-1.094	5.528	-0.198	70.113	0.844
fixed	treatmentmedium	-2.529	8.319	-0.304	70.802	0.762
fixed	treatmenthigh	-0.055	8.113	-0.007	70.113	0.995
fixed	days	0.966	0.107	9.049	112.102	0.000
fixed	treatmentmedium:days	0.473	0.166	2.842	116.319	0.005
fixed	treatmenthigh:days	0.192	0.157	1.225	112.102	0.223

raw area increase by treatment

Fixed Effects: Raw area ~ Treatment x Days

effect	term	estimate	std.error	statistic	df	p.value
fixed	(Intercept)	-0.255	0.077	-3.328	56.632	0.002
fixed	treatmentmedium	-0.095	0.084	-1.132	88.651	0.261
fixed	treatmenthigh	0.022	0.080	0.274	90.189	0.784
fixed	days	0.007	0.001	5.445	111.662	0.000
fixed	initial_area	1.319	0.071	18.455	36.365	0.000
fixed	treatmentmedium:days	0.009	0.002	4.412	116.440	0.000
fixed	treatmenthigh:days	0.000	0.002	0.127	111.662	0.899

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.

Square root transform planar surface area:

# using area_mean
area_sq_lmer <- lmer(sqrt(area_mean) ~ treatment * days + initial_area + (1 | tank/coral), 
                    data = grSumm %>% filter(species == "pstr"))

boundary (singular) fit: see help('isSingular')

# square root raw area: 
tidy(area_sq_lmer, effects = "fixed") |> knitr::kable(digits = 3, caption = "Fixed Effects: Sqrt(Raw area) ~ Treatment x Days") 

Fixed Effects: Sqrt(Raw area) ~ Treatment x Days

effect	term	estimate	std.error	statistic	df	p.value
fixed	(Intercept)	0.401	0.031	13.157	50.789	0.000
fixed	treatmentmedium	-0.029	0.032	-0.908	73.568	0.367
fixed	treatmenthigh	-0.002	0.031	-0.054	74.629	0.957
fixed	days	0.003	0.000	7.934	111.899	0.000
fixed	initial_area	0.590	0.029	20.099	36.881	0.000
fixed	treatmentmedium:days	0.003	0.001	4.154	116.285	0.000
fixed	treatmenthigh:days	0.000	0.001	0.620	111.899	0.537

Square root transformatoin makes treatment effects on PSTR non-significant…

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

term	edf	ref.df	statistic	p.value
s(days):treatmentambient	1.323	1.558	21.591	0.000
s(days):treatmentmedium	2.067	2.390	48.680	0.000
s(days):treatmenthigh	1.000	1.000	28.800	0.000
s(tank)	0.000	3.000	0.000	0.528
s(coral)	24.627	36.000	2.276	0.000

Results. Model fit explains 94% of deviance. All three treatments show significant growth over time, but only medium is significantly higher than ambient.
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

boundary (singular) fit: see help('isSingular')
boundary (singular) fit: see help('isSingular')
boundary (singular) fit: see help('isSingular')

raw area by treatment

Fixed Effects: Raw area ~ Treatment x Days

effect	term	estimate	std.error	statistic	df	p.value
fixed	(Intercept)	-0.258	0.201	-1.282	83.483	0.203
fixed	treatmentmedium	-0.053	0.263	-0.200	102.147	0.842
fixed	treatmenthigh	0.058	0.258	0.226	101.821	0.821
fixed	days	0.009	0.004	2.295	129.486	0.023
fixed	initial_area	1.027	0.012	87.768	40.604	0.000
fixed	treatmentmedium:days	0.014	0.006	2.223	129.837	0.028
fixed	treatmenthigh:days	0.020	0.006	3.370	129.205	0.001

Random Effects: tank and coral

effect	group	term	estimate
ran_pars	coral:tank	sd__(Intercept)	0.462
ran_pars	tank	sd__(Intercept)	0.000
ran_pars	Residual	sd__Observation	0.594

growth rate by treatment

Fixed Effects: Raw area ~ Treatment x Days

effect	term	estimate	std.error	statistic	df	p.value
fixed	(Intercept)	-0.005	0.004	-1.112	98.264	0.269
fixed	treatmentmedium	0.002	0.006	0.263	98.630	0.793
fixed	treatmenthigh	0.011	0.006	1.848	97.998	0.068
fixed	days	0.000	0.000	1.910	81.998	0.060
fixed	treatmentmedium:days	0.000	0.000	-0.180	83.020	0.857
fixed	treatmenthigh:days	0.000	0.000	-1.706	81.263	0.092

Results

All treatments grow significantly over time when considering raw area and days. 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

term	edf	ref.df	statistic	p.value
s(days):treatmentambient	1.007	1.014	5.420	0.021
s(days):treatmentmedium	2.065	2.394	12.238	0.000
s(days):treatmenthigh	1.000	1.000	50.166	0.000
s(tank)	2.354	3.000	724960.042	0.004
s(coral)	39.568	42.000	9730.100	0.245

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

dlabYIIaov <- aov(y ~ treatment, 
                  data = pam_df |>
                    filter(species == "dlab")
                    )
broom::tidy(dlabYIIaov) |> knitr::kable(digits = 3)

term	df	sumsq	meansq	statistic	p.value
treatment	2	0.062	0.031	23.333	0
Residuals	200	0.267	0.001	NA	NA

dlabYIItukey <- glht(dlabYIIaov,
                    linfct = mcp(treatment = "Tukey"))
cld(dlabYIItukey) |> broom::tidy() |> knitr::kable(digits = 3)

treatment	letters
ambient	a
medium	b
high	c

Linear model

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 ' ' 1

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.9601

PSTR

Stats

ANOVA and Tukey

pstrYIIaov <- aov(y ~ treatment, 
                  data = pam_df |>
                    filter(species == "pstr")
                    )
broom::tidy(pstrYIIaov) |> knitr::kable(digits = 3)

term	df	sumsq	meansq	statistic	p.value
treatment	2	0.004	0.002	1.272	0.283
Residuals	189	0.286	0.002	NA	NA

pstrYIItukey <- glht(pstrYIIaov,
                    linfct = mcp(treatment = "Tukey"))
cld(pstrYIItukey) |> broom::tidy() |> knitr::kable(digits = 3)

treatment	letters
ambient	a
medium	a
high	a

#### Linear model

::: {.cell}

```{.r .cell-code  code-fold="true"}
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 ' ' 1

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.1324

:::

Health over time

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, species) |> 
  summarise(n = n(), .groups = "drop") |> 
  group_by(tp, treatment) |> 
  mutate(prop = n / sum(n))

survivor_df |> 
  mutate(health = na_if(health, "")) |> 
  filter(!is.na(health)) |>
  ggplot(aes(x = tp, y = prop, fill = health)) +
  geom_bar(stat = "identity") +
  scale_fill_discrete(labels = c("dead", "healthy", "tissue recession")) +
  facet_grid(species ~ treatment) +
  theme_classic()

At timepoint 0, smilar proportiaons of healthy to tissue recession across treatments and between species, with possible exception of dlab medium treatment. Can see ambient dlab declined in health over time the most. Little change in pstr ambient health. corals in the alkalinity treatments faired better for the most part and seemed to recover from initial tissue recession.