SDS 300 Project Pt. 2

Author

Leah Harries + Lorelei Ing + Destiny Cecchi + Danielle Weiss

Coral Reef % Cover and Carbonate Chemistry

Load Libraries

library(tidyverse)

Warning: package 'tidyverse' was built under R version 4.3.2

Warning: package 'ggplot2' was built under R version 4.3.2

Warning: package 'tidyr' was built under R version 4.3.2

Warning: package 'readr' was built under R version 4.3.2

Warning: package 'purrr' was built under R version 4.3.2

Warning: package 'dplyr' was built under R version 4.3.2

Warning: package 'forcats' was built under R version 4.3.2

Warning: package 'lubridate' was built under R version 4.3.2

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.0
✔ ggplot2   3.4.4     ✔ tibble    3.2.1
✔ lubridate 1.9.3     ✔ tidyr     1.3.0
✔ purrr     1.0.2     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(dplyr)
library(janitor)

Warning: package 'janitor' was built under R version 4.3.3


Attaching package: 'janitor'

The following objects are masked from 'package:stats':

    chisq.test, fisher.test

library(lubridate)
library(leaflet)

Warning: package 'leaflet' was built under R version 4.3.3

library(grDevices)
library(leaflet.extras)

Warning: package 'leaflet.extras' was built under R version 4.3.3

library(broom)

Warning: package 'broom' was built under R version 4.3.2

Data Wrangling

## Read Data
# coral cover data
benthic_cover_florida <- read.csv("C:/Users/Danie/Documents/SDS300_R_WORK_DW/final_proj/benthic_coral_florida.csv")[-1, ]

# carb data
carbonate_chem_atl <- read.csv("C:/Users/Danie/Documents/SDS300_R_WORK_DW/final_proj/carbonate_chem_atlantic.csv")[-1, ]

## Data Cleaning
coral <- clean_names(benthic_cover_florida) |>
  dplyr::select(latitude, longitude, region, mpa_name, year, cover_cat_cd, cover_cat_name, hardbottom_p)

carb <- clean_names(carbonate_chem_atl) |>
  dplyr::select(latitude, longitude, region, year, dic, p_h_measured, aragonite_sat, temperature_c, p_co2) 

# Remove NAs
coral <- coral |>
  filter(!is.na(latitude)) |>
  filter(!is.na(longitude)) 
head(coral)

  latitude longitude region    mpa_name year cover_cat_cd       cover_cat_name
1 24.46499  -81.9823    FLK Unprotected 2016     EUS FAST  Eusmilia fastigiata
2 24.46499  -81.9823    FLK Unprotected 2016     GOR ENCR Encrusting gorgonian
3 24.46499  -81.9823    FLK Unprotected 2016     GOR GORG           Gorgonians
4 24.46499  -81.9823    FLK Unprotected 2016     MEA MEAN Meandrina meandrites
5 24.46499  -81.9823    FLK Unprotected 2016     ORB FAVE  Orbicella faveolata
6 24.46499  -81.9823    FLK Unprotected 2016     SID SIDE  Siderastrea siderea
  hardbottom_p
1            1
2            2
3           13
4            1
5            4
6            4

carb <- carb |>
  filter(!is.na(latitude)) |>
  filter(!is.na(longitude))
head(carb)

     latitude    longitude region year         dic p_h_measured aragonite_sat
1 24.90518333     -80.6266    FLK 2017 2123.014035     7.994548      3.268880
2 24.89826667 -80.61738333    FLK 2017 2140.529652     7.978449      3.187039
3 24.89818333 -80.61606667    FLK 2017 2126.982858     7.986262      3.218789
4     24.8977 -80.61816667    FLK 2017 2145.016273     7.972165      3.162853
5     24.8977 -80.61816667    FLK 2017 2143.278035     7.971244      3.155604
6 24.89663333     -80.6153    FLK 2017 2119.951218     8.014780      3.405095
  temperature_c       p_co2
1          24.7 475.2770518
2          24.7 498.9336095
3          24.7 486.1758318
4          24.8 508.0844039
5          24.8 508.7680168
6          24.7 450.9955019

# Fix Data Types
carb <- carb |>
  mutate(longitude = as.numeric(longitude),
         latitude = as.numeric(latitude),
         temperature_c = as.numeric(temperature_c),
         p_co2 = as.numeric(p_co2),
         dic = as.numeric(dic))
head(coral)

  latitude longitude region    mpa_name year cover_cat_cd       cover_cat_name
1 24.46499  -81.9823    FLK Unprotected 2016     EUS FAST  Eusmilia fastigiata
2 24.46499  -81.9823    FLK Unprotected 2016     GOR ENCR Encrusting gorgonian
3 24.46499  -81.9823    FLK Unprotected 2016     GOR GORG           Gorgonians
4 24.46499  -81.9823    FLK Unprotected 2016     MEA MEAN Meandrina meandrites
5 24.46499  -81.9823    FLK Unprotected 2016     ORB FAVE  Orbicella faveolata
6 24.46499  -81.9823    FLK Unprotected 2016     SID SIDE  Siderastrea siderea
  hardbottom_p
1            1
2            2
3           13
4            1
5            4
6            4

coral <- coral |>
  mutate( longitude = as.numeric(longitude),
          latitude = as.numeric(latitude))
head(carb)

  latitude longitude region year      dic p_h_measured aragonite_sat
1 24.90518 -80.62660    FLK 2017 2123.014     7.994548      3.268880
2 24.89827 -80.61738    FLK 2017 2140.530     7.978449      3.187039
3 24.89818 -80.61607    FLK 2017 2126.983     7.986262      3.218789
4 24.89770 -80.61817    FLK 2017 2145.016     7.972165      3.162853
5 24.89770 -80.61817    FLK 2017 2143.278     7.971244      3.155604
6 24.89663 -80.61530    FLK 2017 2119.951     8.014780      3.405095
  temperature_c    p_co2
1          24.7 475.2771
2          24.7 498.9336
3          24.7 486.1758
4          24.8 508.0844
5          24.8 508.7680
6          24.7 450.9955

Recode Regions, Filter Coral Species, Filter Years

## RECODE REGIONS
# Counted the Regions in coral df
coral |> 
  group_by(region)|>
  summarise(n = n(),
            count = n())

# A tibble: 3 × 3
  region       n count
  <chr>    <int> <int>
1 FLK       3869  3869
2 SEFCRI    3460  3460
3 Tortugas  6272  6272

# Counted the Regions in carb df
carb|> 
  group_by(region)|>
  summarise(n = n(),
            count = n())

# A tibble: 9 × 3
  region     n count
  <chr>  <int> <int>
1 DRTO     115   115
2 DT        43    43
3 FGB       86    86
4 FL       818   818
5 FLK     1404  1404
6 PR        71    71
7 SEFCRI    48    48
8 STX       90    90
9 USVI      90    90

# Recode Region
coral <- coral |>
  mutate(region = recode(region, 'Tortugas' = "DRT")) |>
  filter(region %in% c("FLK", "DRT"))
head(coral)

  latitude longitude region    mpa_name year cover_cat_cd       cover_cat_name
1 24.46499  -81.9823    FLK Unprotected 2016     EUS FAST  Eusmilia fastigiata
2 24.46499  -81.9823    FLK Unprotected 2016     GOR ENCR Encrusting gorgonian
3 24.46499  -81.9823    FLK Unprotected 2016     GOR GORG           Gorgonians
4 24.46499  -81.9823    FLK Unprotected 2016     MEA MEAN Meandrina meandrites
5 24.46499  -81.9823    FLK Unprotected 2016     ORB FAVE  Orbicella faveolata
6 24.46499  -81.9823    FLK Unprotected 2016     SID SIDE  Siderastrea siderea
  hardbottom_p
1            1
2            2
3           13
4            1
5            4
6            4

unique(coral$region)

[1] "FLK" "DRT"

carb <- carb |>
  mutate(region = recode(region, 'FL' = "FLK", 'DRTO' = "DRT", 'DT' = "DRT")) |>
  filter(region %in% c("FLK", "DRT"))
head(carb)

  latitude longitude region year      dic p_h_measured aragonite_sat
1 24.90518 -80.62660    FLK 2017 2123.014     7.994548      3.268880
2 24.89827 -80.61738    FLK 2017 2140.530     7.978449      3.187039
3 24.89818 -80.61607    FLK 2017 2126.983     7.986262      3.218789
4 24.89770 -80.61817    FLK 2017 2145.016     7.972165      3.162853
5 24.89770 -80.61817    FLK 2017 2143.278     7.971244      3.155604
6 24.89663 -80.61530    FLK 2017 2119.951     8.014780      3.405095
  temperature_c    p_co2
1          24.7 475.2771
2          24.7 498.9336
3          24.7 486.1758
4          24.8 508.0844
5          24.8 508.7680
6          24.7 450.9955

unique(carb$region)

[1] "FLK" "DRT"

## FILTER SPECIES
# Count species
count_species <- coral |>
  group_by(cover_cat_cd,cover_cat_name)|>
  summarise(n = n(),
            count = n())|>
  arrange(desc(count))

`summarise()` has grouped output by 'cover_cat_cd'. You can override using the
`.groups` argument.

count_species

# A tibble: 70 × 4
# Groups:   cover_cat_cd [70]
   cover_cat_cd cover_cat_name               n count
   <chr>        <chr>                    <int> <int>
 1 TUR SEDI     Turf Algae with Sediment   789   789
 2 SPO OTHE     Porifera spp               755   755
 3 BAR SUB.     Bare Substrate             729   729
 4 DIC SPE.     Dictyota spp               704   704
 5 GOR GORG     Gorgonians                 698   698
 6 MAC FLES     MacroOtherFleshy           526   526
 7 RHO CRUS     Rhodophyta cru. spp        520   520
 8 GOR ENCR     Encrusting gorgonian       470   470
 9 SID SIDE     Siderastrea siderea        459   459
10 MIL SPE.     Millepora spp              414   414
# ℹ 60 more rows

# Filter out all species more than 20 observations
coral <- coral |>
  group_by(cover_cat_cd,cover_cat_name)|>
  filter(n()>20) |>
  filter(!is.na(cover_cat_name))

# Filter down to only hard corals with > 20 data points (9 hard coral species)
coral <- coral |>
  filter(cover_cat_name %in% c( "Agaricia agaricites", "Colpophyllia natans", "Dichocoenia stokesii",  "Montastraea cavernosa", "Orbicella faveolata", "Porites porites", "Siderastrea siderea", "Stephanocoenia intersepta", "Porites astreoides")) 

head(coral)

# A tibble: 6 × 8
# Groups:   cover_cat_cd, cover_cat_name [5]
  latitude longitude region mpa_name     year cover_cat_cd cover_cat_name       
     <dbl>     <dbl> <chr>  <chr>       <int> <chr>        <chr>                
1     24.5     -82.0 FLK    Unprotected  2016 ORB FAVE     Orbicella faveolata  
2     24.5     -82.0 FLK    Unprotected  2016 SID SIDE     Siderastrea siderea  
3     24.5     -82.0 FLK    Unprotected  2016 STE INTE     Stephanocoenia inter…
4     24.5     -82.0 FLK    Unprotected  2016 DIC STOK     Dichocoenia stokesii 
5     24.5     -81.8 FLK    Unprotected  2016 MON CAVE     Montastraea cavernosa
6     24.5     -81.8 FLK    Unprotected  2016 STE INTE     Stephanocoenia inter…
# ℹ 1 more variable: hardbottom_p <chr>

unique(coral$COVER_CAT_NAME)

Warning: Unknown or uninitialised column: `COVER_CAT_NAME`.

NULL

# Keep only matching years
unique(coral$year)

[1] 2016 2022 2021 2018

unique(carb$year)

[1] 2017 2015 2019 2016 2022 2018 2020 2021

matching_years <- sort(intersect(unique(coral$year), unique(carb$year)))
matching_years

[1] 2016 2018 2021 2022

coral<-  coral[coral$year %in% matching_years, ]
carb <-  carb[carb$year %in% matching_years, ]

# Change levels of years
levels(carb$year) <- c("2016", "2018", "2021", "2022")
levels(coral$year) <- c("2016", "2018", "2021", "2022")

JOIN DATASETS

library(sf)

Warning: package 'sf' was built under R version 4.3.3

Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUE

library(dplyr)

# Convert to sf objects
coral_sf <- st_as_sf(coral, coords = c("longitude", "latitude"), crs = 4326) #n= 1798
carb_sf <- st_as_sf(carb, coords = c("longitude", "latitude"), crs = 4326) #n = 1352

## NEAREST NEIGHBOR

# Extract unique years from the coral dataset
years <- unique(coral$year)

# Make an empty df to store the joined results
joined_df <- data.frame()

# Loop through each year and perform nearest neighbor analysis within that year
for (year in years) {
  # Filter datasets for the current year
  coral_year <- coral_sf[coral_sf$year == year, ]
  carb_year <- carb_sf[carb_sf$year == year, ]
  
  # the nearest neighbor analysis
  nearest_neighbors <- st_nearest_feature(coral_year, carb_year)
  
  # Extract the carbonate chemistry measurements for nearest neighbors
  nearest_carb <- carb_year[nearest_neighbors, ]
  
  # Add the nearest chemistry measurements to the coral cover data
  coral_with_chemistry_year <- cbind(coral_year, nearest_carb)
  
  # Add the result to the final df
  joined_df <- rbind(joined_df, coral_with_chemistry_year)
}

# Reset row numbers
rownames(joined_df) <- NULL

# Update the dataframe
joined_df <- joined_df

#Df check 
unique(joined_df$year)

[1] 2016 2022 2021 2018
Levels: 2016 2018 2021 2022

unique(joined_df$region)

[1] "FLK" "DRT"

unique(joined_df$cover_cat_cd)

[1] "ORB FAVE" "SID SIDE" "STE INTE" "DIC STOK" "MON CAVE" "COL NATA" "POR ASTR"
[8] "AGA AGAR" "POR PORI"

# Now, we can analyze the combined data to understand the relationship between carbonate chemistry and coral cover

Data Visualizations

Map of ΩAr across Florida Keys - in other doc need to combine (with un merged data)

Map of ΩAr, pH, and Temp across Florida Keys - in other doc need to combine (with un merged data)

Statistical Assessments

joined_df

Simple feature collection with 1798 features and 13 fields
Active geometry column: geometry
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -83.10352 ymin: 24.43146 xmax: -80.08719 ymax: 25.74996
Geodetic CRS:  WGS 84
First 10 features:
   region    mpa_name year cover_cat_cd            cover_cat_name hardbottom_p
1     FLK Unprotected 2016     ORB FAVE       Orbicella faveolata            4
2     FLK Unprotected 2016     SID SIDE       Siderastrea siderea            4
3     FLK Unprotected 2016     STE INTE Stephanocoenia intersepta            1
4     FLK Unprotected 2016     DIC STOK      Dichocoenia stokesii            1
5     FLK Unprotected 2016     MON CAVE     Montastraea cavernosa            1
6     FLK Unprotected 2016     STE INTE Stephanocoenia intersepta            1
7     FLK  West_Sambo 2016     MON CAVE     Montastraea cavernosa            1
8     FLK  West_Sambo 2016     SID SIDE       Siderastrea siderea            1
9     FLK  West_Sambo 2016     STE INTE Stephanocoenia intersepta            1
10    FLK  West_Sambo 2016     COL NATA       Colpophyllia natans            1
   region.1 year.1    dic p_h_measured aragonite_sat temperature_c p_co2
1       FLK   2016 2018.0        8.003          3.88        30.492 448.5
2       FLK   2016 2018.0        8.003          3.88        30.492 448.5
3       FLK   2016 2018.0        8.003          3.88        30.492 448.5
4       FLK   2016 2018.0        8.003          3.88        30.492 448.5
5       FLK   2016 2039.6        8.006          3.85        29.634 447.6
6       FLK   2016 2039.6        8.006          3.85        29.634 447.6
7       FLK   2016 2053.9        8.035          3.65        26.287 416.3
8       FLK   2016 2053.9        8.035          3.65        26.287 416.3
9       FLK   2016 2053.9        8.035          3.65        26.287 416.3
10      FLK   2016 2053.3        8.059          3.74        25.728 392.1
                     geometry                 geometry.1
1   POINT (-81.9823 24.46499) POINT (-81.82954 24.45697)
2   POINT (-81.9823 24.46499) POINT (-81.82954 24.45697)
3   POINT (-81.9823 24.46499) POINT (-81.82954 24.45697)
4   POINT (-81.9823 24.46499) POINT (-81.82954 24.45697)
5   POINT (-81.7548 24.46302) POINT (-81.71873 24.47505)
6   POINT (-81.7548 24.46302) POINT (-81.71873 24.47505)
7  POINT (-81.72699 24.47505) POINT (-81.71912 24.47733)
8  POINT (-81.72699 24.47505) POINT (-81.71912 24.47733)
9  POINT (-81.72699 24.47505) POINT (-81.71912 24.47733)
10 POINT (-81.70926 24.50635) POINT (-81.71482 24.47673)

# additive multiple regression

# How does the percentage of hard coral cover (hardbottom_p) vary over time (year) in each region (region)?
coral_lm1 <- lm(hardbottom_p ~ year + region, data = joined_df)
summary(coral_lm1)


Call:
lm(formula = hardbottom_p ~ year + region, data = joined_df)

Residuals:
   Min     1Q Median     3Q    Max 
-3.097 -1.442 -1.135  0.558 37.558 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.25599    0.15612  14.450  < 2e-16 ***
year2018     0.08407    0.19747   0.426    0.670    
year2021     0.18556    0.24626   0.753    0.451    
year2022    -0.12133    0.19342  -0.627    0.531    
regionFLK    0.84059    0.15930   5.277 1.47e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.087 on 1793 degrees of freedom
Multiple R-squared:  0.01557,   Adjusted R-squared:  0.01338 
F-statistic: 7.091 on 4 and 1793 DF,  p-value: 1.158e-05

#            Pr(>|t|)    
# (Intercept)  < 2e-16 ***
# year2018       0.670    
# year2021       0.451    
# year2022       0.531    
# regionFLK   1.47e-07 ***
# F-statistic: 7.091 on 4 and 1793 DF,  p-value: 1.158e-05

# How does the percentage of hard coral cover (hardbottom_p) vary over time (year) in each region (region) - for each species?
coral_lm2 <- lm(hardbottom_p ~ year + region + cover_cat_name, data = joined_df)
summary(coral_lm2)


Call:
lm(formula = hardbottom_p ~ year + region + cover_cat_name, data = joined_df)

Residuals:
   Min     1Q Median     3Q    Max 
-5.368 -1.490 -0.517  0.510 36.841 

Coefficients:
                                        Estimate Std. Error t value Pr(>|t|)
(Intercept)                              1.13249    0.30613   3.699 0.000223
year2018                                 0.15562    0.19035   0.818 0.413735
year2021                                 0.15281    0.23693   0.645 0.519022
year2022                                 0.10136    0.18888   0.537 0.591563
regionFLK                                0.90368    0.15506   5.828 6.64e-09
cover_cat_nameColpophyllia natans        1.63062    0.44434   3.670 0.000250
cover_cat_nameDichocoenia stokesii       1.87377    0.53476   3.504 0.000470
cover_cat_nameMontastraea cavernosa      1.42356    0.32322   4.404 1.12e-05
cover_cat_nameOrbicella faveolata        3.23044    0.35597   9.075  < 2e-16
cover_cat_namePorites astreoides         0.35197    0.31254   1.126 0.260250
cover_cat_namePorites porites            0.06277    0.38792   0.162 0.871480
cover_cat_nameSiderastrea siderea        1.23613    0.29804   4.148 3.52e-05
cover_cat_nameStephanocoenia intersepta  0.22878    0.32426   0.706 0.480568
                                           
(Intercept)                             ***
year2018                                   
year2021                                   
year2022                                   
regionFLK                               ***
cover_cat_nameColpophyllia natans       ***
cover_cat_nameDichocoenia stokesii      ***
cover_cat_nameMontastraea cavernosa     ***
cover_cat_nameOrbicella faveolata       ***
cover_cat_namePorites astreoides           
cover_cat_namePorites porites              
cover_cat_nameSiderastrea siderea       ***
cover_cat_nameStephanocoenia intersepta    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.964 on 1785 degrees of freedom
Multiple R-squared:  0.09602,   Adjusted R-squared:  0.08995 
F-statistic:  15.8 on 12 and 1785 DF,  p-value: < 2.2e-16

# F-statistic:  15.8 on 12 and 1785 DF,  p-value: < 2.2e-16

# add in chemistry stats
# use augment() for additive models in ggplot
# additive multiple regression
# 1 categorical

# How does carbonate chemistry effect hard coral cover over time in each region? 
# exploring models:
coral_lm3 <- lm(hardbottom_p ~ year + region + p_h_measured, data = joined_df)
summary(coral_lm3)


Call:
lm(formula = hardbottom_p ~ year + region + p_h_measured, data = joined_df)

Residuals:
   Min     1Q Median     3Q    Max 
-3.240 -1.453 -1.091  0.547 37.547 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -42.59455   24.17125  -1.762  0.07821 .  
year2018       0.08330    0.19734   0.422  0.67298    
year2021       0.01751    0.26223   0.067  0.94678    
year2022      -0.10558    0.19347  -0.546  0.58534    
regionFLK      0.67643    0.18212   3.714  0.00021 ***
p_h_measured   5.59928    3.01755   1.856  0.06368 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.085 on 1792 degrees of freedom
Multiple R-squared:  0.01746,   Adjusted R-squared:  0.01472 
F-statistic: 6.369 on 5 and 1792 DF,  p-value: 7.15e-06

# F-statistic: 6.369 on 5 and 1792 DF,  p-value: 7.15e-06

coral_lm4 <- lm(hardbottom_p ~ year + region + p_h_measured + dic, data = joined_df)
summary(coral_lm4)


Call:
lm(formula = hardbottom_p ~ year + region + p_h_measured + dic, 
    data = joined_df)

Residuals:
   Min     1Q Median     3Q    Max 
-3.239 -1.453 -1.096  0.547 37.547 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -4.193e+01  2.549e+01  -1.645 0.100161    
year2018      8.125e-02  1.990e-01   0.408 0.683026    
year2021      1.719e-02  2.623e-01   0.066 0.947766    
year2022     -1.004e-01  2.036e-01  -0.493 0.622184    
regionFLK     6.718e-01  1.908e-01   3.520 0.000442 ***
p_h_measured  5.573e+00  3.035e+00   1.837 0.066432 .  
dic          -2.218e-04  2.695e-03  -0.082 0.934403    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.085 on 1791 degrees of freedom
Multiple R-squared:  0.01747,   Adjusted R-squared:  0.01417 
F-statistic: 5.306 on 6 and 1791 DF,  p-value: 1.954e-05

# F-statistic: 5.306 on 6 and 1791 DF,  p-value: 1.954e-05

coral_lm5 <- lm(hardbottom_p ~ year + region + p_h_measured + dic + aragonite_sat, data = joined_df)
summary(coral_lm5)


Call:
lm(formula = hardbottom_p ~ year + region + p_h_measured + dic + 
    aragonite_sat, data = joined_df)

Residuals:
   Min     1Q Median     3Q    Max 
-3.889 -1.462 -1.042  0.492 37.538 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)   
(Intercept)   -38.925191  25.437766  -1.530  0.12614   
year2018        0.455018   0.228794   1.989  0.04688 * 
year2021       -0.166689   0.267553  -0.623  0.53335   
year2022       -0.011399   0.204885  -0.056  0.95564   
regionFLK       0.541843   0.194371   2.788  0.00536 **
p_h_measured    7.096764   3.061805   2.318  0.02057 * 
dic            -0.004763   0.003023  -1.576  0.11530   
aragonite_sat  -1.601079   0.488008  -3.281  0.00105 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.077 on 1790 degrees of freedom
Multiple R-squared:  0.02334,   Adjusted R-squared:  0.01952 
F-statistic: 6.111 on 7 and 1790 DF,  p-value: 4.493e-07

# F-statistic: 6.111 on 7 and 1790 DF,  p-value: 4.493e-07

coral_lm6 <- lm(hardbottom_p ~ year + region + p_h_measured + dic + aragonite_sat + temperature_c, data = joined_df)
summary(coral_lm6)


Call:
lm(formula = hardbottom_p ~ year + region + p_h_measured + dic + 
    aragonite_sat + temperature_c, data = joined_df)

Residuals:
   Min     1Q Median     3Q    Max 
-3.779 -1.459 -0.937  0.489 37.541 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)   107.604752  46.787741   2.300 0.021571 *  
year2018        0.576103   0.230280   2.502 0.012447 *  
year2021       -0.703242   0.303003  -2.321 0.020404 *  
year2022        0.625715   0.266303   2.350 0.018900 *  
regionFLK       0.580323   0.193950   2.992 0.002808 ** 
p_h_measured  -10.043812   5.520064  -1.820 0.069001 .  
dic            -0.008143   0.003146  -2.589 0.009714 ** 
aragonite_sat   0.188357   0.683453   0.276 0.782890    
temperature_c  -0.317490   0.085211  -3.726 0.000201 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.066 on 1789 degrees of freedom
Multiple R-squared:  0.03086,   Adjusted R-squared:  0.02652 
F-statistic:  7.12 on 8 and 1789 DF,  p-value: 2.614e-09

# F-statistic:  7.12 on 8 and 1789 DF,  p-value: 2.614e-09

coral_lm7 <- lm(hardbottom_p ~ year + region + p_h_measured + dic + aragonite_sat + temperature_c + p_co2, data = joined_df)
summary(coral_lm7)


Call:
lm(formula = hardbottom_p ~ year + region + p_h_measured + dic + 
    aragonite_sat + temperature_c + p_co2, data = joined_df)

Residuals:
   Min     1Q Median     3Q    Max 
-4.083 -1.478 -1.010  0.482 37.542 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)   
(Intercept)   134.301573  49.788143   2.697  0.00705 **
year2018        0.514694   0.233514   2.204  0.02764 * 
year2021       -0.712905   0.302943  -2.353  0.01872 * 
year2022        0.587656   0.267306   2.198  0.02805 * 
regionFLK       0.597809   0.194194   3.078  0.00211 **
p_h_measured  -11.934646   5.648801  -2.113  0.03476 * 
dic            -0.018696   0.007446  -2.511  0.01213 * 
aragonite_sat   3.879706   2.457716   1.579  0.11461   
temperature_c  -0.817895   0.331182  -2.470  0.01362 * 
p_co2           0.023682   0.015146   1.564  0.11810   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.065 on 1788 degrees of freedom
Multiple R-squared:  0.03218,   Adjusted R-squared:  0.02731 
F-statistic: 6.606 on 9 and 1788 DF,  p-value: 2.491e-09

# F-statistic: 6.606 on 9 and 1788 DF,  p-value: 2.491e-09

# ANOVAS - still working on 
anova1 <- aov(aragonite_sat~region + year, data = joined_df)
summary(anova1)

              Df Sum Sq Mean Sq F value Pr(>F)    
region         1   0.00   0.000       0  0.989    
year           3  28.99   9.662     333 <2e-16 ***
Residuals   1793  52.02   0.029                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(anova1)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = aragonite_sat ~ region + year, data = joined_df)

$region
                diff         lwr        upr     p adj
FLK-DRT -0.000107943 -0.01609473 0.01587884 0.9894357

$year
                 diff         lwr         upr     p adj
2018-2016  0.25791625  0.23001890  0.28581360 0.0000000
2021-2016 -0.07840671 -0.11231450 -0.04449893 0.0000000
2022-2016 -0.01162758 -0.03873244  0.01547727 0.6876635
2021-2018 -0.33632296 -0.37016084 -0.30248508 0.0000000
2022-2018 -0.26954383 -0.29656119 -0.24252648 0.0000000
2022-2021  0.06677913  0.03359158  0.09996668 0.0000015

anova2 <- aov(p_co2 ~ cover_cat_cd, data= joined_df)
summary(anova2)

               Df  Sum Sq Mean Sq F value   Pr(>F)    
cover_cat_cd    8   96996   12124   5.043 3.25e-06 ***
Residuals    1789 4300969    2404                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(anova2)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = p_co2 ~ cover_cat_cd, data = joined_df)

$cover_cat_cd
                         diff        lwr        upr     p adj
COL NATA-AGA AGAR -23.0494605 -45.716743 -0.3821783 0.0427397
DIC STOK-AGA AGAR   0.4167899 -26.934516 27.7680960 1.0000000
MON CAVE-AGA AGAR  -8.2991467 -24.712898  8.1146042 0.8210280
ORB FAVE-AGA AGAR -13.9395674 -32.112487  4.2333521 0.2943523
POR ASTR-AGA AGAR  -0.4780264 -16.483599 15.5275462 1.0000000
POR PORI-AGA AGAR  13.4983660  -6.383544 33.3802763 0.4677354
SID SIDE-AGA AGAR  -9.0101940 -24.277339  6.2569507 0.6601961
STE INTE-AGA AGAR  -9.9999756 -26.559463  6.5595121 0.6309554
DIC STOK-COL NATA  23.4662503  -6.479676 53.4121772 0.2662478
MON CAVE-COL NATA  14.7503137  -5.696588 35.1972159 0.3798113
ORB FAVE-COL NATA   9.1098931 -12.774338 30.9941246 0.9336301
POR ASTR-COL NATA  22.5714340   2.450725 42.6921433 0.0148676
POR PORI-COL NATA  36.5478265  13.224985 59.8706682 0.0000433
SID SIDE-COL NATA  14.0392665  -5.499165 33.5776977 0.3854930
STE INTE-COL NATA  13.0494849  -7.514591 33.6135608 0.5641785
MON CAVE-DIC STOK  -8.7159366 -34.257355 16.8254822 0.9796482
ORB FAVE-DIC STOK -14.3563573 -41.062309 12.3495944 0.7652576
POR ASTR-DIC STOK  -0.8948163 -26.175861 24.3862281 1.0000000
POR PORI-DIC STOK  13.0815761 -14.815435 40.9785869 0.8753934
SID SIDE-DIC STOK  -9.4269839 -34.247108 15.3931403 0.9608181
STE INTE-DIC STOK -10.4167655 -36.052083 15.2185517 0.9420419
ORB FAVE-MON CAVE  -5.6404207 -20.954634  9.6737925 0.9674036
POR ASTR-MON CAVE   7.8211203  -4.845706 20.4879462 0.6015944
POR PORI-MON CAVE  21.7975127   4.489697 39.1053284 0.0030620
SID SIDE-MON CAVE  -0.7110473 -12.430938 11.0088432 1.0000000
STE INTE-MON CAVE  -1.7008289 -15.060722 11.6590647 0.9999837
POR ASTR-ORB FAVE  13.4615410  -1.414354 28.3374363 0.1129586
POR PORI-ORB FAVE  27.4379334   8.453617 46.4222497 0.0002619
SID SIDE-ORB FAVE   4.9293734  -9.148965 19.0077120 0.9760952
STE INTE-ORB FAVE   3.9395918 -11.530720 19.4099034 0.9971264
POR PORI-POR ASTR  13.9763925  -2.944826 30.8976104 0.2022136
SID SIDE-POR ASTR  -8.5321676 -19.673216  2.6088813 0.2964867
STE INTE-POR ASTR  -9.5219491 -22.377061  3.3331624 0.3423124
SID SIDE-POR PORI -22.5085600 -38.733078 -6.2840424 0.0005859
STE INTE-POR PORI -23.4983416 -40.944427 -6.0522562 0.0010062
STE INTE-SID SIDE  -0.9897816 -12.912921 10.9333575 0.9999994

anova3 <- aov(dic + temperature_c ~ region, data=joined_df)
summary(anova3)

              Df  Sum Sq Mean Sq F value Pr(>F)    
region         1  111866  111866   125.2 <2e-16 ***
Residuals   1796 1604784     894                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(anova3)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = dic + temperature_c ~ region, data = joined_df)

$region
             diff       lwr       upr p adj
FLK-DRT -16.00594 -18.81157 -13.20032     0