Environmental and local habitat variables as predictors of trophic interactions in subtidal rocky reefs along the SE Pacific coast

Clear the environment
Load packages
Part 1: Importing and examining data
- 1a) Import data
- 1b) Inspect structure of data frames
Part 2: Manipulating datasets
Part 3: Latitudinal patterns in explanatory variables
- 3a): Explore latitudinal relationships
Part 4: Latitudinal patterns in consumption
- 4a): Predation vs latitude
- 4b): Herbivory vs latitude
Part 5: Predation intensity: Effect of site-specific variables
Part 6: Herbivory intensity: Effect of site specific variables
- 6a) Model selection
Part 7: Visualisation of best fitted predictors
Part 9: Supplementary tables
- 9a) Size estimations summary table
Part 10: R Session Information

Clear the environment

rm(list = ls())

Load packages

Before loading the packages install them using install()

library(tidyverse)
library(summarytools)
library(dplyr)
library(readr)
library(vegan)
library(patchwork)
library(pander)
library(ggeffects)
library(performance)
library(knitr)
library(car)
library(corrplot)
library(MASS)
library(MuMIn)
library(geosphere)
library(mgcv)

#If Issues installing packages
options("install.lock"=FALSE)

Part 1: Importing and examining data

In this part the data frames were inspected in order to check the variable names and types. In some situations the data set was slightly modified for further manipulation of them.

1a) Import data

Environmental1 <- read.csv("Data/Environmental_variables.csv")
Metadata <- read.csv("Data/Metadata_lat2018-2019.csv")
Herbivory <- read.csv("Data/Herbivory_latitudinal2018-2019.csv")
Benthic <- read.csv("Data/Percent_cover.csv")
Phlorotannins <- read.csv("Data/Phlorotannins.csv")
Predation <- read.csv("Data/Predation_latitudinal2018-2019.csv")
Size <- read.csv("Data/Size_estimations.csv")
Transects <- read.csv("Data/Transect2018-2019.csv")

1b) Inspect structure of data frames

Site metadata

view(dfSummary(Metadata))

#Rename temperature variables to avoid confusion with replicate measures in Environmental data
Metadata <- rename(Metadata, Temperature_site = Temperature_C)

Environmental data

view(dfSummary(Environmental1))

Environmental <- Environmental1 %>%
  rename(D_oxygen_percent = X., D_oxygen_mgL = D_oxygen) %>% #Rename oxygen percent variable
  dplyr::select(-Date,-Campaign,-Latitude) #Remove variables that are in Metadata

Benthic cover

view(dfSummary(Benthic))
Benthic <- Benthic %>%
  dplyr::select(-Date) #Remove variables that are in Metadata

Herbivory

view(dfSummary(Herbivory))
Herbivory <- Herbivory %>%
  dplyr::select(-Date,-Campaign,-Season) %>%# Remove variables that are in Metadata
  mutate(Consumed_proportion = (Initial_length_mm-Final_Length_mm)/Initial_length_mm)

It was necessary to create major categories for Bait

Herbivory <- mutate(Herbivory, Bait2 = case_when(Bait == "Lessonia_spicata" ~ "Lessonia sp.",
                                                 Bait == "Lessonia_berteroana" ~ "Lessonia sp."))

Predation

view(dfSummary(Predation))
Predation <- Predation %>%
  dplyr::select(-Date,-Campaign,-Season)  %>%# Remove variables that are in Metadata
mutate(Consumed_proportion = (Initial_Individuals-Remaining_individuals)/Initial_Individuals, Consumed_individuals = Initial_Individuals-Remaining_individuals) #Add Consumed_proportion column

Create major categories for Bait (Squidpops vs Crabs) and Assay time (2 hours and 24 hours)

Predation <- mutate(Predation, Bait2 = case_when(Bait == "Squidpops" ~ "Squidpops",
                                                 Bait != "Squidpops" ~ "Crabs"), Assay_time2 = case_when(Assay_time == "24" ~ "24",
                                                 Bait != "24" ~ "2"))

Phlorotannins

view(dfSummary(Phlorotannins))
Phlorotannins <- Phlorotannins %>%
  dplyr::select(-Date,-Latitude) # Remove variables that are in metadata

Size

view(dfSummary(Size))
Size <- Size %>%
  dplyr::select(-Date,-Campaign) # Remove variables that are in Metadata

Transects

view(dfSummary(Transects))
Transects <- Transects %>%
  dplyr::select(-Date,-Campaign) # Remove variables that are in Metadata

Part 2: Manipulating datasets

In this section data sets are manipulated for analyses. New data sets with means and combined with metadata and consumption are created.

2a) Environmental data

Join Environmental data with metadata

Environmental_metadata <- Environmental %>%
    left_join(Metadata, by = "Site")

Join Environmental data with consumption

First summarise environmental data per site

Environmental_2 <- Environmental_metadata %>%
  group_by(Site, Latitude) %>%
  dplyr::select(-D_oxygen_percent, -Conductivity_Ms, -C_f, -C_i, -H_Visibility) #Keep only variables that will be used in the analysis

Environmental_2$Day_length <- daylength(Environmental_2$Latitude, Environmental_2$Date)


Environmental_sum_long <- Environmental_2 %>%
  dplyr::select(Site, Latitude, Temperature, D_oxygen_mgL, Salinity_ppt, V_Visibility, Day_length) %>%
  gather(Env_Variable, Value, Temperature:Day_length) %>%
  group_by(Site, Latitude, Env_Variable)%>%
  summarise(Mean_value = mean(Value, na.rm=T),n = n(), sd = sd(Value,na.rm=T), se = sd/sqrt(n)) 

Environmental_sum <- Environmental_sum_long %>%
  dplyr::select(Site, Env_Variable, Mean_value) %>%
pivot_wider(names_from = Env_Variable, values_from = Mean_value,
              values_fill = 0)

Then, add environmental data to consumption
Predation vs environmental data sets

Predation_env <- Predation %>%
    left_join(Environmental_sum, by = "Site")

Herbivory vs environmental data sets

Herbivory_env <- Herbivory %>%
  left_join(Environmental_sum, by = "Site")

2b): Benthic data

Join benthic data with metadata

Benthic_metadata <- Benthic  %>%
    left_join(Metadata, by = "Site")

Benthic_metadata_wide <- Benthic_metadata %>%
  group_by(Site, Latitude, Type) %>%
  pivot_wider(names_from = Type, values_from = Coverage,
              values_fill = 0)

Calculate mean, sd, and se by group

Benthic_mean <- Benthic_metadata %>%
  group_by(Site, Latitude, Type)%>%
  summarise(Mean_cover = mean(Coverage),n = n(), sd = sd(Coverage), se = sd/sqrt(n))

Join benthic cover data with consumption

First, summarise benthic cover per site and change it to the wide format

Benthic_sum_wide <- Benthic %>%
  group_by(Site, Type) %>%
  summarise(cover_mean = mean(Coverage)) %>%
  pivot_wider(names_from = Type, values_from = cover_mean,
              values_fill = 0)
Benthic_sum_wide_var <- dplyr::select(Benthic_sum_wide, Site, Kelp, Bare_rock)

Then, add benthic data to consumption
Predation_env + Benthic summarised data

Predation_2 <- Predation_env %>%
    left_join(Benthic_sum_wide_var, by = "Site")

Herbivory_env + Benthic summarised data

Herbivory_2 <- Herbivory_env %>%
  left_join(Benthic_sum_wide_var, by = "Site")

2c): Fish length

Join Size data with metadata

Size_metadata <- Size %>%
left_join(Metadata, by = c("Site"))

Calculate mean, sd and se per group

Size_mean_sp<- Size_metadata %>%
  group_by(Site, Latitude, Species, Type)%>%
  summarise(mean_Size_sp = mean(Size_cm))
Size_mean <- Size_mean_sp %>%
   group_by(Site, Latitude, Type)%>%
  summarise(mean_Size = mean(mean_Size_sp),n = n(), sd = sd(mean_Size_sp), se = sd/sqrt(n))

Calculate mean of fish size and add it to consumption

Size_mean_wide <- Size_mean_sp %>%
  filter(Type == "Fish") %>%
  group_by(Site, Type)%>%
  summarise(mean_Size = mean(mean_Size_sp)) %>%
  pivot_wider(names_from = Type, values_from = mean_Size,
              values_fill = 0) %>%
  rename(Fish_size = Fish)

Add fish length data to consumption
Predation_env + fish size data

Predation_3 <- Predation_2 %>%
    left_join(Size_mean_wide, by = "Site")

2d) Species list (abundances and richness)

Add metadata to Transects data set

#Join transects with metadata
Transect_metadata <- Transects %>%
left_join(Metadata, by = c("Site"))

#Summary of species from transects
Transect_sum <- Transect_metadata %>%
  filter(Method == "Transect") %>% #Select only data taken with visual transect method
  group_by(Site, Latitude, Census, Species) %>%
  summarise(Total_abundance = sum(Individuals))

Calculate abundances

Calculate abundances by species

First, change data frame to the wide format

Transect_wide <- Transect_sum %>%
   ungroup()  %>%
  pivot_wider(names_from = Species, values_from = Total_abundance,
              values_fill = 0)

Convert back to long format with Species and Abundance as columns. This will add rows with zeroes for all the species missing in each video

Transect_long <- Transect_wide %>%
  gather(Species, Abundance, Anisotremus_scapularis:Labrisomus_philippii)

#Obtain mean abundance per site and method
Transect_mean <- Transect_long %>%
      group_by(Site, Latitude, Species) %>%
     summarise(Mean_abundance = mean(Abundance),n = n(), sd = sd(Abundance), se = sd/sqrt(n))

Calculate abundances by group

Follow the same steps as abundances by type

TransectGroup_sum <- Transect_metadata %>%
  group_by(Site, Latitude, Census, Group) %>%
  summarise(Total_abundance = sum(Individuals))

TransectGroup_wide <- TransectGroup_sum %>%
  ungroup()  %>%
  pivot_wider(names_from = Group, values_from = Total_abundance,
              values_fill = 0) 
TransectGroup_wide <- TransectGroup_wide%>%
  mutate(Herbivores = Herbivore_fish + Herbivore_gastropod + Urchin + Omnivorous_crab, Invertebrate_predators = Predatory_gastropod+ Snail + Predatory_crab + Sea_star + Omnivorous_crab)


TransectGroup_long <- TransectGroup_wide %>%
  gather(Group, Abundance, Herbivore_fish:Invertebrate_predators)

TransectGroup_mean <- TransectGroup_long %>%
  group_by(Site, Latitude, Group) %>%
  summarise(Total_abundance = mean(Abundance),n = n(), sd = sd(Abundance), se = sd/sqrt(n))

Calculate Species Richness

First, add a column of presence when abundance is higher than 0.

Transect_long <- Transect_long %>%
  mutate(Presence = if_else(Abundance > 0,1,0))

Calculate sum of all species in the same site and census.

Richness_transect <- Transect_long %>%
  group_by(Site, Latitude, Census) %>%
  summarise(Richness = sum(Presence))

#Calculate mean by site  
Richness_mean <- Richness_transect %>%
  group_by(Site, Latitude)%>%
  summarise(Mean_SR = mean(Richness),n = n(), sd = sd(Richness), se = sd/sqrt(n))

Add abundances to consumption data

Calculate means and change it to the wide format

#Per Group
Group_abundances_mean <- TransectGroup_long%>%
  filter(Group == "Invertebrate_predators" | Group == "Predatory_fish"| Group == "Herbivores")%>%
  group_by(Site, Group)%>%
  summarise(Mean_Abundance = mean(Abundance)) %>%
  pivot_wider(names_from = Group, values_from = Mean_Abundance,
              values_fill = 0)

Then, add abundances data to consumption

Predation_env + group abundance summarised data

#Select only the needed columns
Group_abundances_mean_pred <- dplyr::select(Group_abundances_mean, Site, Invertebrate_predators, Predatory_fish)

Predation_4 <- Predation_3 %>%
    left_join(Group_abundances_mean_pred, by = "Site")

Herbivory_env + Benthic summarised data

#Select only the needed columns
Group_abundances_mean_herb <- dplyr::select(Group_abundances_mean, Site, Herbivores)

Herbivory_3 <- Herbivory_2 %>%
  left_join(Group_abundances_mean_herb, by = "Site")

Add species richness to consumption data

Calculate means and change it to the wide format

Richness_meanpred <- Richness_transect %>%
  group_by(Site)%>%
  summarise(Mean_SR = mean(Richness))

Then, add abundances data to consumption

Predation_4 + species richness summarised data

Predation_5 <- Predation_4 %>%
    left_join(Richness_meanpred, by = "Site")

Herbivory_3 + species richness summarised data

Herbivory_4 <- Herbivory_3 %>%
  left_join(Richness_meanpred, by = "Site")

2e): Community composition

Summarise abundances of species per site and standardise

Transect_sum <- Transect_long %>%
  group_by(Site, Latitude, Species) %>%
  summarise(Total_abundance = mean(Abundance))

## `summarise()` has grouped output by 'Site', 'Latitude'. You can override using
## the `.groups` argument.

Transect_wide2 <- Transect_sum %>%
ungroup() %>%
pivot_wider(names_from = Species, values_from = Total_abundance, values_fill = 0)

Transect_val <- Transect_wide2 %>%
  dplyr::select(- Site, -Latitude) 

Transect_stand <- decostand(Transect_val, method = "total")

Calculate Bray-Curtis dissimilariry

Transect_comp <- vegdist(Transect_stand, "bray")

Run Principal Coordinate Analysis (PcoA) to obtain the component explaining the highers variation

pcoa <- cmdscale(Transect_comp, eig = TRUE, add = TRUE)
pcoa$eig

##  [1]  1.740530e+00  1.329554e+00  8.438524e-01  6.203963e-01  5.740214e-01
##  [6]  3.517972e-01  1.671264e-01  7.777519e-02  5.878784e-02  3.878297e-02
## [11] -7.537429e-17 -2.615278e-16

(pcoa$eig/sum(pcoa$eig))*100 #proportion of variation

##  [1]  2.999557e+01  2.291298e+01  1.454260e+01  1.069165e+01  9.892446e+00
##  [6]  6.062727e+00  2.880187e+00  1.340345e+00  1.013125e+00  6.683695e-01
## [11] -1.298969e-15 -4.507061e-15

#convert pcoa results into data frame that can be plotted
pcoa_df <- data.frame(pcoa$points)
colnames(pcoa_df) <- c("PCo1", "PCo2")
pcoa_df <- pcoa_df %>%
  dplyr::select(-c(PCo2))
pcoa_df$Site <- factor(Transect_wide2$Site)
pcoa_df$Latitude <- Transect_wide2$Latitude

PCo 1 represents 30% of the variation community composition

Add PCos to consumption data

Predation_4 + PCo

pcoa_df1<- dplyr::select(pcoa_df, -Latitude)
Predation_6 <- Predation_5 %>%
    left_join(pcoa_df1, by = "Site")

Herbivory_3 + PCo

Herbivory_5 <- Herbivory_4  %>%
    left_join(pcoa_df1, by = "Site")

2f): Phlorotannins

Add metadata

Phlorotannins_metadata <- Phlorotannins %>%
left_join(Metadata, by = c("Site"))

Summarise phlorotannins per site

Phlorotannins_mean <- Phlorotannins_metadata %>%
  group_by(Site, Latitude) %>%
  summarise(phloro_mean = mean(Phlorotannins), n = n(), sd = sd(Phlorotannins), se = sd/sqrt(n))

Herbivory + phlorotannins

Phlorotannins_sum <- Phlorotannins %>%
  group_by(Site) %>%
  summarise(phloro_mean = mean(Phlorotannins))

Herbivory_6 <- Herbivory_5  %>%
    left_join(Phlorotannins_sum, by = "Site")

Part 3: Latitudinal patterns in explanatory variables

Latitudinal patterns were analysed for the recorded explanatory variables using linear models.

3a): Explore latitudinal relationships

Set general theme for plots

Theme1 <- theme_bw() + 
    theme(panel.grid.major = element_line(colour = "white"),
          panel.grid.minor = element_line(colour = "white"), axis.text.y = element_text(size = 10), axis.title.y = element_text(size = 9), axis.text.x = element_text(size = 10), axis.title.x =     element_text(size = 10))

For these analyses we used the approach and code of model selection for environmental and spatial gradients provided by Anderson et al. (2022)

1) Temperature vs latitude

a. Analyses

For future plotting

Environmental_2$Latitude <- Environmental_2$Latitude*-1
#change sign of latitude

# First, examine a simple scatter plot (always a good idea!)
   plot(Environmental_2$Latitude, Environmental_2$Temperature, xlab = "Latitude", ylab = "Temperature", las=1)
   mygrey <- grey(level = 0.65, alpha = 0.4)
   points(Environmental_2$Latitude, Environmental_2$Temperature, pch = 21, col = mygrey, bg = mygrey)

T_mod1 <- mgcv::gam(Temperature ~ s(Latitude), data= Environmental_2, family = gaussian())
T_mod2 <- mgcv::gam(Temperature ~ s(Latitude), data= Environmental_2, family = Gamma())
T_mod3 <- mgcv::gam(Temperature ~ s(Latitude), data= Environmental_2, family = inverse.gaussian())

AICc(T_mod1,T_mod2,T_mod3)

##              df     AICc
## T_mod1 10.49323 92.95936
## T_mod2 10.79361 87.83173
## T_mod3 10.89346 95.94008

T_mod2 is the best fitted model

summary(T_mod2)

## 
## Family: Gamma 
## Link function: inverse 
## 
## Formula:
## Temperature ~ s(Latitude)
## 
## Parametric coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.0669789  0.0003234   207.1   <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(Latitude) 8.794  8.988 107.4  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.945   Deviance explained = 96.7%
## GCV = 0.0013495  Scale est. = 0.0010856  n = 48

b. Plot predicted relationship

new_data <- data.frame(Latitude = seq(min(Environmental_2$Latitude), max(Environmental_2$Latitude), length.out = 100))
predictions <- predict(T_mod2, newdata = new_data, type = "link", se.fit = TRUE)

# Combine predicted values with new_data
ilink <- family(T_mod2)$linkinv
pred <- cbind(new_data, predictions)
pred <- transform(pred, lwr_ci = ilink(fit - (2 * se.fit)),
                        upr_ci = ilink(fit + (2 * se.fit)),
                        fitted = ilink(fit))

# Plot the GAM model predictions
Temperature_plot <- ggplot() +
  geom_jitter(data = Environmental_2, aes(x= Latitude, y=Temperature), width = 0.08, height = 0.03, size = 1, alpha = 0.6) +
  geom_ribbon(data = pred, aes(x = Latitude, y = fitted, ymin = lwr_ci, ymax = upr_ci), alpha = 0.2) +
  geom_line(data = pred, aes(x = Latitude, y = fitted), linewidth = 0.6) + 
  Theme1 +
labs( y="Temperature (°C)", x="Latitude (°S)") +
theme(axis.title.x=element_blank(),
        axis.text.x=element_blank())

par(mfrow = c(2, 2))
gam.check(T_mod2)

## 
## Method: GCV   Optimizer: outer newton
## full convergence after 7 iterations.
## Gradient range [2.08852e-11,2.08852e-11]
## (score 0.001349529 & scale 0.001085617).
## Hessian positive definite, eigenvalue range [1.289977e-05,1.289977e-05].
## Model rank =  10 / 10 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##               k'  edf k-index p-value  
## s(Latitude) 9.00 8.79    0.76   0.055 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

2) Salinity vs latitude

a. Analyses

# First, examine a simple scatter plot (always a good idea!)
   plot(Environmental_2$Latitude, Environmental_2$Salinity_ppt, xlab = "Latitude", ylab = "Salinity", las=1)
   mygrey <- grey(level = 0.65, alpha = 0.4)
   points(Environmental_2$Latitude, Environmental_2$Salinity_ppt, pch = 21, col = mygrey, bg = mygrey)

S_mod1 <- gam(Salinity_ppt ~ s(Latitude), data= Environmental_2, family = gaussian())
S_mod2 <- gam(Salinity_ppt ~ s(Latitude), data= Environmental_2, family = Gamma())
S_mod3 <- gam(Salinity_ppt ~ s(Latitude), data= Environmental_2, family = inverse.gaussian())

AICc(S_mod1, S_mod3, S_mod2)

##              df     AICc
## S_mod1 8.835490 18.62663
## S_mod3 8.896160 27.17283
## S_mod2 8.875632 19.53924

summary(S_mod1)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Salinity_ppt ~ s(Latitude)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 33.73125    0.03678   917.2   <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(Latitude) 6.835  7.953 8.194   2e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.569   Deviance explained = 63.2%
## GCV = 0.077583  Scale est. = 0.064918  n = 48

b. Plot predicted relationship

new_data <- data.frame(Latitude = seq(min(Environmental_2$Latitude), max(Environmental_2$Latitude), length.out = 100))
predictions <- predict(S_mod1, newdata = new_data, type = "link", se.fit = TRUE)

# Combine predicted values with new_data
ilink <- family(S_mod1)$linkinv
pred <- cbind(new_data, predictions)
pred <- transform(pred, lwr_ci = ilink(fit - (2 * se.fit)),
                        upr_ci = ilink(fit + (2 * se.fit)),
                        fitted = ilink(fit))

# Plot the GAM model predictions
Salinity_plot <- ggplot() +
  geom_jitter(data = Environmental_2, aes(x= Latitude, y=Salinity_ppt), width = 0.08, height = 0.05, size = 1, alpha = 0.6) +
  geom_ribbon(data = pred, aes(x = Latitude, y = fitted, ymin = lwr_ci, ymax = upr_ci), alpha = 0.2) +
  geom_line(data = pred, aes(x = Latitude, y = fitted), linewidth = 0.6) + 
  Theme1 +
labs( y="Salinity (ppt)", x="Latitude (°S)") +
theme(axis.title.x=element_blank(),
        axis.text.x=element_blank())

par(mfrow = c(2, 2))
gam.check(S_mod1)

## 
## Method: GCV   Optimizer: magic
## Smoothing parameter selection converged after 9 iterations.
## The RMS GCV score gradient at convergence was 1.068683e-07 .
## The Hessian was positive definite.
## Model rank =  10 / 10 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##               k'  edf k-index p-value
## s(Latitude) 9.00 6.84    1.37    0.99

3) Oxygen vs latitude

a. Analyses

# First, examine a simple scatter plot (always a good idea!)
   plot(Environmental_2$Latitude, Environmental_2$D_oxygen_mgL, xlab = "Latitude", ylab = "DO", las=1)
   mygrey <- grey(level = 0.65, alpha = 0.4)
   points(Environmental_2$Latitude, Environmental_2$D_oxygen_mgL, pch = 21, col = mygrey, bg = mygrey)

Ox_mod1 <- gam(D_oxygen_mgL ~ s(Latitude), data= Environmental_2, family = gaussian())
Ox_mod2 <- gam(D_oxygen_mgL ~ s(Latitude), data= Environmental_2, family = Gamma())
Ox_mod3 <- gam(D_oxygen_mgL ~ s(Latitude), data= Environmental_2, family = inverse.gaussian())

AICc(Ox_mod1, Ox_mod2, Ox_mod3)

##               df     AICc
## Ox_mod1 8.700722 226.6735
## Ox_mod2 9.661964 230.3983
## Ox_mod3 4.230413 210.4661

Two options may be possible based on AICc. Very low df in Ox_mod 3

summary( Ox_mod1)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## D_oxygen_mgL ~ s(Latitude)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    6.381      0.322   19.82   <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(Latitude) 6.701  7.838 2.711  0.0183 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.279   Deviance explained = 38.1%
## GCV = 5.9291  Scale est. = 4.9779    n = 48

Deviance explained in Ox_mod3_12 is too low compared to the other. Then, Ox_mod1_12 was chosen

b. Plot predicted relationship

new_data <- data.frame(Latitude = seq(min(Environmental_2$Latitude), max(Environmental_2$Latitude), length.out = 100))
predictions <- predict(Ox_mod1, newdata = new_data, type = "link", se.fit = TRUE)

# Combine predicted values with new_data
ilink <- family(Ox_mod1)$linkinv
pred <- cbind(new_data, predictions)
pred <- transform(pred, lwr_ci = ilink(fit - (2 * se.fit)),
                        upr_ci = ilink(fit + (2 * se.fit)),
                        fitted = ilink(fit))

# Plot the GAM model predictions
Oxygen_plot <- ggplot() +
  geom_jitter(data = Environmental_2, aes(x= Latitude, y=D_oxygen_mgL), width = 0.08, height = 0.07, size = 1, alpha = 0.6) +
  geom_ribbon(data = pred, aes(x = Latitude, y = fitted, ymin = lwr_ci, ymax = upr_ci), alpha = 0.2) +
  geom_line(data = pred, aes(x = Latitude, y = fitted), linewidth = 0.6) + 
  Theme1 +
labs( y="Dissolved oxygen (mgL)", x="Latitude (°S)") +
theme(axis.title.x=element_text(size = 10),
        axis.text.x=element_text(size = 10), axis.text.y = element_text(size = 10), axis.title.y = element_text(size = 8))

par(mfrow = c(2, 2))
gam.check(Ox_mod1)

## 
## Method: GCV   Optimizer: magic
## Smoothing parameter selection converged after 7 iterations.
## The RMS GCV score gradient at convergence was 3.970196e-05 .
## The Hessian was positive definite.
## Model rank =  10 / 10 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##              k' edf k-index p-value  
## s(Latitude) 9.0 6.7    0.73    0.04 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

4) Visibility vs latitude

a. Analyses

# First, examine a simple scatter plot (always a good idea!)
   plot(Environmental_2$Latitude, Environmental_2$V_Visibility, xlab = "Latitude", ylab = "Visibility", las=1)
   mygrey <- grey(level = 0.65, alpha = 0.4)
   points(Environmental_2$Latitude, Environmental_2$V_Visibility, pch = 21, col = mygrey, bg = mygrey)

V_mod1 <- gam(V_Visibility ~ s(Latitude),gamma=1.4, data= Environmental_2, family = gaussian())
V_mod2 <- gam(V_Visibility ~ s(Latitude), data= Environmental_2, family = inverse.gaussian())
V_mod3 <- gam(V_Visibility ~ s(Latitude), data= Environmental_2, family = Gamma())

AICc(V_mod1, V_mod2, V_mod3)

##              df     AICc
## V_mod1 8.590546 187.0415
## V_mod2 7.856608 195.2612
## V_mod3 8.317530 191.4637

V_mod1_12 is the best fitted

summary(V_mod1)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## V_Visibility ~ s(Latitude)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   5.8542     0.2136   27.41   <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(Latitude) 6.591  7.741 5.224 0.000246 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =   0.45   Deviance explained = 52.7%
## GCV = 3.0404  Scale est. = 2.1894    n = 48

b. Plot predicted relationship

new_data <- data.frame(Latitude = seq(min(Environmental_2$Latitude), max(Environmental_2$Latitude), length.out = 100))
predictions <- predict(V_mod1, newdata = new_data, type = "link", se.fit = TRUE)

# Combine predicted values with new_data
ilink <- family(V_mod1)$linkinv
pred <- cbind(new_data, predictions)
pred <- transform(pred, lwr_ci = ilink(fit - (2 * se.fit)),
                        upr_ci = ilink(fit + (2 * se.fit)),
                        fitted = ilink(fit))

# Plot the GAM model predictions
Visibility_plot <- ggplot() +
  geom_jitter(data = Environmental_2, aes(x= Latitude, y=V_Visibility), width = 0.08, height = 0.08, size = 1, alpha = 0.6) +
  geom_ribbon(data = pred, aes(x = Latitude, y = fitted, ymin = lwr_ci, ymax = upr_ci), alpha = 0.2) +
  geom_line(data = pred, aes(x = Latitude, y = fitted), linewidth = 0.6) + 
  Theme1 +
labs( y="Visibility (m)", x="Latitude (°S)") +
theme(axis.title.x=element_blank(),
        axis.text.x=element_blank())

par(mfrow = c(2, 2))
gam.check(V_mod1)

## 
## Method: GCV   Optimizer: magic
## Smoothing parameter selection converged after 9 iterations.
## The RMS GCV score gradient at convergence was 4.595236e-06 .
## The Hessian was positive definite.
## Model rank =  10 / 10 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##               k'  edf k-index p-value   
## s(Latitude) 9.00 6.59    0.68    0.01 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

5) PCOA vs latitude

a. Analyses

pcoa_df$Latitude <- pcoa_df$Latitude*-1

# First, examine a simple scatter plot (always a good idea!)
   plot(pcoa_df$Latitude, pcoa_df$PCo1, xlab = "Latitude", ylab = "PCoA1", las=1)
   mygrey <- grey(level = 0.65, alpha = 0.4)
   points(pcoa_df$Latitude, pcoa_df$PCo1, pch = 21, col = mygrey, bg = mygrey)

PCO_mod1 <- gam(PCo1~ s(Latitude), gamma= 1.4, data= pcoa_df,  family = gaussian())

AICc(PCO_mod1)

## [1] 15.22301

summary(PCO_mod1)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## PCo1 ~ s(Latitude)
## 
## Parametric coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.205e-17  9.808e-02       0        1
## 
## Approximate significance of smooth terms:
##              edf Ref.df     F p-value  
## s(Latitude) 1.09  1.173 4.554  0.0634 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.271   Deviance explained = 34.3%
## GCV = 0.1667  Scale est. = 0.11543   n = 12

par(mfrow = c(2, 2))
gam.check(PCO_mod1)

## 
## Method: GCV   Optimizer: magic
## Smoothing parameter selection converged after 9 iterations.
## The RMS GCV score gradient at convergence was 3.007551e-06 .
## The Hessian was positive definite.
## Model rank =  10 / 10 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##               k'  edf k-index p-value
## s(Latitude) 9.00 1.09    1.29    0.72

No significant effect

# Plot the GAM model predictions
PCOA_plot <- ggplot() +
  geom_jitter(data = pcoa_df, aes(x= Latitude, y=PCo1), width = 0, height = 0, size = 1, alpha = 0.6) +
  Theme1 +
labs( y="PCoA1", x="Latitude (°S)")

6) Kelp vs latitude

a. Analyses

For further plotting

Benthic_metadata_wide$Latitude <- Benthic_metadata_wide$Latitude*-1

# First, examine a simple scatter plot (always a good idea!)
   plot(Benthic_metadata_wide$Latitude, Benthic_metadata_wide$Kelp, xlab = "Latitude", ylab = "Kelp", las=1)
   mygrey <- grey(level = 0.65, alpha = 0.4)
   points(Benthic_metadata_wide$Latitude, Benthic_metadata_wide$Kelp, pch = 21, col = mygrey, bg = mygrey)

Kelp_mod1 <- mgcv::gam(Kelp ~ s(Latitude), gamma = 1.4, data = Benthic_metadata_wide, family = gaussian())
Kelp_mod2 <- mgcv::gam(Kelp ~ s(Latitude), gamma =1.4, data = Benthic_metadata_wide, family = tw())

AICc(Kelp_mod2, Kelp_mod1)

##                 df     AICc
## Kelp_mod2 4.640437 476.8968
## Kelp_mod1 4.217000 663.0220

summary(Kelp_mod2)

## 
## Family: Tweedie(p=1.365) 
## Link function: log 
## 
## Formula:
## Kelp ~ s(Latitude)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    2.351      0.156   15.07   <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(Latitude) 1.367   1.64 13.1 4.69e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =   0.33   Deviance explained = 24.4%
## -REML = 169.53  Scale est. = 7.6607    n = 78

par(mfrow = c(2, 2))
gam.check(Kelp_mod2)

## 
## Method: REML   Optimizer: outer newton
## full convergence after 8 iterations.
## Gradient range [-2.196392e-05,1.485091e-05]
## (score 169.5299 & scale 7.660684).
## Hessian positive definite, eigenvalue range [0.06262817,51.51969].
## Model rank =  10 / 10 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##               k'  edf k-index p-value
## s(Latitude) 9.00 1.37    0.79    0.12

b. Plot predicted relationship

new_data <- data.frame(Latitude = seq(min(Benthic_metadata_wide$Latitude), max(Benthic_metadata_wide$Latitude), length.out = 100))
predictions <- predict(Kelp_mod2, newdata = new_data, type = "link", se.fit = TRUE)

# Combine predicted values with new_data
ilink <- family(Kelp_mod2)$linkinv
pred <- cbind(new_data, predictions)
pred <- transform(pred, lwr_ci = ilink(fit - (2 * se.fit)),
                        upr_ci = ilink(fit + (2 * se.fit)),
                        fitted = ilink(fit))

# Plot the GAM model predictions
Kelp_plot <- ggplot() +
  geom_jitter(data = Benthic_metadata_wide, aes(x= Latitude, y=Kelp), width = 0.08, height = 0.11, size = 1, alpha = 0.6) +
  geom_ribbon(data = pred, aes(x = Latitude, y = fitted, ymin = lwr_ci, ymax = upr_ci), alpha = 0.2) +
  geom_line(data = pred, aes(x = Latitude, y = fitted), linewidth = 0.6) + 
  Theme1 +
labs( y="Kelp cover (%)", x="Latitude (°S)") +
theme(axis.title.x=element_blank(),
        axis.text.x=element_blank())

7) Bare_rock vs latitude

a. Analyses

# First, examine a simple scatter plot (always a good idea!)
   plot(Benthic_metadata_wide$Latitude, Benthic_metadata_wide$Bare_rock, xlab = "Latitude", ylab = "Bare_rock", las=1)
   mygrey <- grey(level = 0.65, alpha = 0.4)
   points(Benthic_metadata_wide$Latitude, Benthic_metadata_wide$Bare_rock, pch = 21, col = mygrey, bg = mygrey)

BR_mod1 <- mgcv::gam(Bare_rock ~ s(Latitude), gamma= 1.4, data = Benthic_metadata_wide, family = gaussian())
BR_mod2 <- mgcv::gam(Bare_rock ~ s(Latitude), gamma=1.4, data = Benthic_metadata_wide, family = tw())

AICc(BR_mod1, BR_mod2)

##                df     AICc
## BR_mod1  9.728495 706.3690
## BR_mod2 11.606470 637.8918

summary(BR_mod2)

## 
## Family: Tweedie(p=1.3) 
## Link function: log 
## 
## Formula:
## Bare_rock ~ s(Latitude)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.20356    0.08913   35.94   <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(Latitude) 7.926  8.606 12.08  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.598   Deviance explained = 62.8%
## -REML = 238.88  Scale est. = 4.7114    n = 78

par(mfrow = c(2, 2))
gam.check(BR_mod2)

## 
## Method: REML   Optimizer: outer newton
## full convergence after 7 iterations.
## Gradient range [-1.240802e-07,7.992647e-07]
## (score 238.878 & scale 4.711392).
## Hessian positive definite, eigenvalue range [1.777301,45.39948].
## Model rank =  10 / 10 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##               k'  edf k-index p-value  
## s(Latitude) 9.00 7.93    0.75   0.015 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

b. Plot predicted relationship

new_data <- data.frame(Latitude = seq(min(Benthic_metadata_wide$Latitude), max(Benthic_metadata_wide$Latitude), length.out = 100))
predictions <- predict(BR_mod2, newdata = new_data, type = "link", se.fit = TRUE)

# Combine predicted values with new_data
ilink <- family(BR_mod2)$linkinv
pred <- cbind(new_data, predictions)
pred <- transform(pred, lwr_ci = ilink(fit - (2 * se.fit)),
                        upr_ci = ilink(fit + (2 * se.fit)),
                        fitted = ilink(fit))

# Plot the GAM model predictions
Bare_plot <- ggplot() +
  geom_jitter(data = Benthic_metadata_wide, aes(x= Latitude, y=Bare_rock), width = 0.08, height = 0.11, size = 1, alpha = 0.6) +
  geom_ribbon(data = pred, aes(x = Latitude, y = fitted, ymin = lwr_ci, ymax = upr_ci), alpha = 0.2) +
  geom_line(data = pred, aes(x = Latitude, y = fitted), linewidth = 0.6) + 
  Theme1 +
labs( y="Bare rock cover(%)", x="Latitude (°S)") +
theme(axis.title.x=element_blank(),
        axis.text.x=element_blank())

8) Predatory fish abundance vs latitude

TransectGroup_wide$Latitude <- TransectGroup_wide$Latitude*-1

# First, examine a simple scatter plot (always a good idea!)
   plot(TransectGroup_wide$Latitude, TransectGroup_wide$Predatory_fish, xlab = "Latitude", ylab = "PF", las=1)
   mygrey <- grey(level = 0.65, alpha = 0.4)
   points(TransectGroup_wide$Latitude, TransectGroup_wide$Predatory_fish, pch = 21, col = mygrey, bg = mygrey)

PF_mod1 <- mgcv::gam(Predatory_fish ~ s(Latitude), gamma=1.4, data = TransectGroup_wide, family =poisson)
PF_mod2 <- mgcv::gam(Predatory_fish ~ s(Latitude), gamma=1.4, data = TransectGroup_wide, family =quasipoisson)
PF_mod3 <- mgcv::gam(Predatory_fish ~ s(Latitude), gamma=1.4, data = TransectGroup_wide, family =nb)

AICc(PF_mod2, PF_mod1, PF_mod3)

##               df      AICc
## PF_mod2 4.624668        NA
## PF_mod1 9.516203 1174.2311
## PF_mod3 6.155803  252.0494

summary(PF_mod3)

## 
## Family: Negative Binomial(1.038) 
## Link function: log 
## 
## Formula:
## Predatory_fish ~ s(Latitude)
## 
## Parametric coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   3.4929     0.2019    17.3   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##               edf Ref.df Chi.sq  p-value    
## s(Latitude) 3.351  4.156  30.09 6.78e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.0722   Deviance explained = 49.6%
## -REML = 89.213  Scale est. = 1         n = 26

par(mfrow = c(2, 2))
gam.check(PF_mod3)

## 
## Method: REML   Optimizer: outer newton
## full convergence after 2 iterations.
## Gradient range [-1.29827e-07,5.38212e-07]
## (score 89.21338 & scale 1).
## Hessian positive definite, eigenvalue range [0.604836,7.703663].
## Model rank =  10 / 10 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##               k'  edf k-index p-value
## s(Latitude) 9.00 3.35    1.03    0.67

b. Plot predicted relationship

new_data <- data.frame(Latitude = seq(min(TransectGroup_wide$Latitude), max(TransectGroup_wide$Latitude), length.out = 100))
predictions <- predict(PF_mod3, newdata = new_data, type = "link", se.fit = TRUE)

# Combine predicted values with new_data
ilink <- family(PF_mod3)$linkinv
pred <- cbind(new_data, predictions)
pred <- transform(pred, lwr_ci = ilink(fit - (2 * se.fit)),
                        upr_ci = ilink(fit + (2 * se.fit)),
                        fitted = ilink(fit))

# Plot the GAM model predictions
PredFish_plot <- ggplot() +
  geom_jitter(data = TransectGroup_wide, aes(x= Latitude, y=Predatory_fish), width = 0.12, height = 0.5, size = 1, alpha = 0.6) +
  geom_ribbon(data = pred, aes(x = Latitude, y = fitted, ymin = lwr_ci, ymax = upr_ci), alpha = 0.2) +
  geom_line(data = pred, aes(x = Latitude, y = fitted), linewidth = 0.6) + 
  Theme1 +
labs( y="Predatory Fish", x="Latitude (°S)") +
theme(axis.title.x=element_blank(),
        axis.text.x=element_blank())

9) Herbivorous fish abundance vs latitude

a. Analyses

# First, examine a simple scatter plot (always a good idea!)
   plot(TransectGroup_wide$Latitude, TransectGroup_wide$Herbivores, xlab = "Latitude", ylab = "Herbivores", las=1)
   mygrey <- grey(level = 0.65, alpha = 0.4)
   points(TransectGroup_wide$Latitude, TransectGroup_wide$Herbivores, pch = 21, col = mygrey, bg = mygrey)

Herb_mod1 <- mgcv::gam(Herbivores ~ s(Latitude), gamma=1.4, data = TransectGroup_wide, family =poisson())
Herb_mod2 <- mgcv::gam(Herbivores ~ s(Latitude), gamma=1.4, data = TransectGroup_wide, family =quasipoisson())
Herb_mod3 <- mgcv::gam(Herbivores ~ s(Latitude), gamma=1.4, data = TransectGroup_wide, family =nb())

AICc(Herb_mod1, Herb_mod2, Herb_mod3)

##                 df      AICc
## Herb_mod1 9.983415 4296.9233
## Herb_mod2 8.545300        NA
## Herb_mod3 7.826680  333.2693

summary(Herb_mod3)

## 
## Family: Negative Binomial(0.454) 
## Link function: log 
## 
## Formula:
## Herbivores ~ s(Latitude)
## 
## Parametric coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   4.9115     0.3178   15.45   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##               edf Ref.df Chi.sq  p-value    
## s(Latitude) 4.757  5.827  30.12 3.74e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.195   Deviance explained =   50%
## -REML = 117.44  Scale est. = 1         n = 26

par(mfrow = c(2, 2))
gam.check(Herb_mod3)

## 
## Method: REML   Optimizer: outer newton
## full convergence after 6 iterations.
## Gradient range [-3.297643e-09,5.452203e-09]
## (score 117.4409 & scale 1).
## Hessian positive definite, eigenvalue range [0.9303036,8.671097].
## Model rank =  10 / 10 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##               k'  edf k-index p-value
## s(Latitude) 9.00 4.76    0.79    0.24

b. Plot predicted relationship

new_data <- data.frame(Latitude = seq(min(TransectGroup_wide$Latitude), max(TransectGroup_wide$Latitude), length.out = 100))
predictions <- predict(Herb_mod3, newdata = new_data, type = "link", se.fit = TRUE)

# Combine predicted values with new_data
ilink <- family(Herb_mod3)$linkinv
pred <- cbind(new_data, predictions)
pred <- transform(pred, lwr_ci = ilink(fit - (2 * se.fit)),
                        upr_ci = ilink(fit + (2 * se.fit)),
                        fitted = ilink(fit))

# Plot the GAM model predictions
Herb_plot <- ggplot() +
  geom_jitter(data = TransectGroup_wide, aes(x= Latitude, y=Herbivores), width = 0.2, height = 0.5, size = 1, alpha = 0.6) +
  geom_ribbon(data = pred, aes(x = Latitude, y = fitted, ymin = lwr_ci, ymax = upr_ci), alpha = 0.2) +
  geom_line(data = pred, aes(x = Latitude, y = fitted), linewidth = 0.6) + 
  Theme1 +
labs( y="Herbivores", x="Latitude (°S)")

10) Invertebrate predators abundance vs latitude

# First, examine a simple scatter plot (always a good idea!)
   plot(TransectGroup_wide$Latitude, TransectGroup_wide$Invertebrate_predators, xlab = "Latitude", ylab = "Invertebrate pred", las=1)
   mygrey <- grey(level = 0.65, alpha = 0.4)
   points(TransectGroup_wide$Latitude, TransectGroup_wide$Invertebrate_predators, pch = 21, col = mygrey, bg = mygrey)

Inv_mod1 <- mgcv::gam(Invertebrate_predators ~ s(Latitude), gamma=1.4, data = TransectGroup_wide, family =poisson())
Inv_mod2 <- mgcv::gam(Invertebrate_predators ~ s(Latitude), gamma=1.4, data = TransectGroup_wide, family =quasipoisson())
Inv_mod3 <- mgcv::gam(Invertebrate_predators ~ s(Latitude), gamma=1.4, data = TransectGroup_wide, family = nb())

AICc(Inv_mod2, Inv_mod1, Inv_mod3)

##                df     AICc
## Inv_mod2 7.690064       NA
## Inv_mod1 9.885648 660.0753
## Inv_mod3 8.229799 272.5043

summary(Inv_mod3)

## 
## Family: Negative Binomial(1.273) 
## Link function: log 
## 
## Formula:
## Invertebrate_predators ~ s(Latitude)
## 
## Parametric coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   3.7830     0.1779   21.26   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##               edf Ref.df Chi.sq p-value    
## s(Latitude) 5.119   6.23  35.49  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.473   Deviance explained = 66.2%
## -REML = 97.611  Scale est. = 1         n = 26

par(mfrow = c(2, 2))
gam.check(Inv_mod3)

## 
## Method: REML   Optimizer: outer newton
## full convergence after 4 iterations.
## Gradient range [-5.190965e-05,4.311985e-05]
## (score 97.61055 & scale 1).
## Hessian positive definite, eigenvalue range [0.9545727,6.48156].
## Model rank =  10 / 10 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##               k'  edf k-index p-value
## s(Latitude) 9.00 5.12    0.86    0.17

b. Plot predicted relationship

new_data <- data.frame(Latitude = seq(min(TransectGroup_wide$Latitude), max(TransectGroup_wide$Latitude), length.out = 100))
predictions <- predict(Inv_mod3, newdata = new_data, type = "link", se.fit = TRUE)

# Combine predicted values with new_data
ilink <- family(Inv_mod3)$linkinv
pred <- cbind(new_data, predictions)
pred <- transform(pred, lwr_ci = ilink(fit - (2 * se.fit)),
                        upr_ci = ilink(fit + (2 * se.fit)),
                        fitted = ilink(fit))

# Plot the GAM model predictions
Inv_plot <- ggplot() +
  geom_jitter(data = TransectGroup_wide, aes(x= Latitude, y=Invertebrate_predators), width = 0.08, height = 0.5, size = 1, alpha = 0.6) +
  geom_ribbon(data = pred, aes(x = Latitude, y = fitted, ymin = lwr_ci, ymax = upr_ci), alpha = 0.2) +
  geom_line(data = pred, aes(x = Latitude, y = fitted), linewidth = 0.6) + 
  Theme1 +
labs( y="Invertebrate predators", x="Latitude (°S)") +
theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(), axis.title.y = element_text(size = 8))

11) Species richness vs latitude

Richness_transect$Latitude <- Richness_transect$Latitude*-1

# First, examine a simple scatter plot (always a good idea!)
   plot(Richness_transect$Latitude, Richness_transect$Richness, xlab = "Latitude", ylab = "Richness", las=1)
   mygrey <- grey(level = 0.65, alpha = 0.4)
   points(Richness_transect$Latitude, Richness_transect$Richness, pch = 21, col = mygrey, bg = mygrey)

R_mod1 <- mgcv::gam(Richness ~ s(Latitude), gamma=1.4, data = Richness_transect, family =poisson())
R_mod2 <- mgcv::gam(Richness ~ s(Latitude), gamma=1.4,data = Richness_transect, family =nb())

AICc(R_mod2, R_mod1)

##              df     AICc
## R_mod2 3.013601 136.7409
## R_mod1 4.934598 130.7241

summary(R_mod1)

## 
## Family: poisson 
## Link function: log 
## 
## Formula:
## Richness ~ s(Latitude)
## 
## Parametric coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  2.07428    0.07004   29.62   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##               edf Ref.df Chi.sq p-value  
## s(Latitude) 3.935  4.855  8.917   0.099 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.302   Deviance explained = 43.2%
## UBRE = 0.19305  Scale est. = 1         n = 26

par(mfrow = c(2, 2))
gam.check(R_mod1)

## 
## Method: UBRE   Optimizer: outer newton
## full convergence after 4 iterations.
## Gradient range [-6.400339e-12,-6.400339e-12]
## (score 0.1930469 & scale 1).
## Hessian positive definite, eigenvalue range [0.06283876,0.06283876].
## Model rank =  10 / 10 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##               k'  edf k-index p-value
## s(Latitude) 9.00 3.93    0.98    0.36

No significant effect

# Plot the GAM model predictions
Richness_plot <- ggplot() +
  geom_jitter(data = Richness_transect, aes(x= Latitude, y=Richness), width = 0.08, height = 0.5, size = 1, alpha = 0.6) +
  Theme1 +
labs( y="Species Richness", x="Latitude (°S)") +
theme(axis.title.y = element_text(size = 8))

12) Fish length vs latitude

Size_fish <- filter(Size_mean_sp, Type == "Fish")
Size_fish$Latitude <- Size_fish$Latitude*-1

# First, examine a simple scatter plot (always a good idea!)
   plot(Size_fish$Latitude, Size_fish$mean_Size_sp, xlab = "Latitude", ylab = "Fish Length", las=1)
   mygrey <- grey(level = 0.65, alpha = 0.4)
   points(Size_fish$Latitude, Size_fish$mean_Size_sp, pch = 21, col = mygrey, bg = mygrey)

Length_mod1 <- mgcv::gam(mean_Size_sp ~ s(Latitude), gamma=1.4, data = Size_fish, family =gaussian())
Length_mod2 <- mgcv::gam(mean_Size_sp ~ s(Latitude), gamma=1.4, data = Size_fish, family =Gamma())
Length_mod3 <- mgcv::gam(mean_Size_sp ~ s(Latitude), gamma=1.4, data = Size_fish, family =inverse.gaussian())

AICc(Length_mod1, Length_mod2, Length_mod3)

##                   df     AICc
## Length_mod1 3.716434 876.0274
## Length_mod2 3.823501 861.9641
## Length_mod3 3.801599 851.5481

summary(Length_mod3)

## 
## Family: inverse.gaussian 
## Link function: 1/mu^2 
## 
## Formula:
## mean_Size_sp ~ s(Latitude)
## 
## Parametric coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.0020493  0.0001609   12.74   <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(Latitude) 1.802  2.244 3.408  0.0305 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.0578   Deviance explained = 7.13%
## GCV = 0.010315  Scale est. = 0.0082538  n = 119

par(mfrow = c(2, 2))
gam.check(Length_mod3)

## 
## Method: GCV   Optimizer: outer newton
## full convergence after 4 iterations.
## Gradient range [-7.194329e-14,-7.194329e-14]
## (score 0.01031498 & scale 0.008253778).
## Hessian positive definite, eigenvalue range [9.707093e-05,9.707093e-05].
## Model rank =  10 / 10 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##              k' edf k-index p-value
## s(Latitude) 9.0 1.8    1.12    0.94

new_data <- data.frame(Latitude = seq(min(Size_fish$Latitude), max(Size_fish$Latitude), length.out = 100))
predictions <- predict(Length_mod3, newdata = new_data, type = "link", se.fit = TRUE)

# Combine predicted values with new_data
ilink <- family(Length_mod3)$linkinv
pred <- cbind(new_data, predictions)
pred <- transform(pred, lwr_ci = ilink(fit - (2 * se.fit)),
                        upr_ci = ilink(fit + (2 * se.fit)),
                        fitted = ilink(fit))

# Plot the GAM model predictions
Size_plot <- ggplot() +
  geom_jitter(data = Size_fish, aes(x= Latitude, y=mean_Size_sp), width = 0.08, height = 0.5, size = 1, alpha = 0.6) +
  geom_ribbon(data = pred, aes(x = Latitude, y = fitted, ymin = lwr_ci, ymax = upr_ci), alpha = 0.2) +
  geom_line(data = pred, aes(x = Latitude, y = fitted), linewidth = 0.6) + 
  Theme1 +
labs( y="Fish length (cm)", x="Latitude (°S)") +
theme(axis.title.x=element_blank(),
        axis.text.x=element_blank())

13) Phlorotannins vs latitude

Phlorotannins_metadata$Latitude <- Phlorotannins_metadata$Latitude*-1

# First, examine a simple scatter plot (always a good idea!)
   plot(Phlorotannins_metadata$Latitude, Phlorotannins_metadata$Phlorotannins, xlab = "Latitude", ylab = "phlorotannins", las=1)
   mygrey <- grey(level = 0.65, alpha = 0.4)
   points(Phlorotannins_metadata$Latitude, Phlorotannins_metadata$Phlorotannins, pch = 21, col = mygrey, bg = mygrey)

Phloro_mod1 <- mgcv::gam(Phlorotannins ~ s(Latitude), gamma=1.4, data = Phlorotannins_metadata, family =gaussian())
Phloro_mod2 <- mgcv::gam(Phlorotannins ~ s(Latitude), gamma=1.4, data = Phlorotannins_metadata, family =Gamma())
Phloro_mod3 <- mgcv::gam(Phlorotannins ~ s(Latitude), gamma=1.4,  data = Phlorotannins_metadata, family =inverse.gaussian())

AICc(Phloro_mod1, Phloro_mod2, Phloro_mod3)

##                    df     AICc
## Phloro_mod1  9.380438 934.2137
## Phloro_mod2  8.836817 817.7763
## Phloro_mod3 10.190366 821.3613

summary(Phloro_mod2)

## 
## Family: Gamma 
## Link function: inverse 
## 
## Formula:
## Phlorotannins ~ s(Latitude)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.078330   0.003997    19.6   <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(Latitude) 6.837  7.933 15.99  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.452   Deviance explained = 56.3%
## GCV = 0.32923  Scale est. = 0.2625    n = 120

par(mfrow = c(2, 2))
gam.check(Phloro_mod2)

## 
## Method: GCV   Optimizer: outer newton
## full convergence after 4 iterations.
## Gradient range [3.958715e-07,3.958715e-07]
## (score 0.3292323 & scale 0.2624964).
## Hessian positive definite, eigenvalue range [0.003156813,0.003156813].
## Model rank =  10 / 10 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##               k'  edf k-index p-value
## s(Latitude) 9.00 6.84    0.86    0.12

new_data <- data.frame(Latitude = seq(min(Phlorotannins_metadata$Latitude), max(Phlorotannins_metadata$Latitude), length.out = 100))
predictions <- predict(Phloro_mod2, newdata = new_data, type = "link", se.fit = TRUE)

# Combine predicted values with new_data
ilink <- family(Phloro_mod2)$linkinv
pred <- cbind(new_data, predictions)
pred <- transform(pred, lwr_ci = ilink(fit - (2 * se.fit)),
                        upr_ci = ilink(fit + (2 * se.fit)),
                        fitted = ilink(fit))

# Plot the GAM model predictions
Phloro_plot <- ggplot() +
  geom_jitter(data = Phlorotannins_metadata, aes(x= Latitude, y=Phlorotannins), width = 0.08, height = 0.5, size = 1, alpha = 0.6) +
  geom_ribbon(data = pred, aes(x = Latitude, y = fitted, ymin = lwr_ci, ymax = upr_ci), alpha = 0.2) +
  geom_line(data = pred, aes(x = Latitude, y = fitted), linewidth = 0.6) + 
  Theme1 +
labs( y="Phlorotannins (mg"~g^-1*" DW)", x="Latitude (°S)")

Final plots

Figure 2 - Abiotic predictors

Temperature_plot + Salinity_plot + Oxygen_plot + plot_layout(ncol=1)

Figure 2: Latitudinal variation of abiotic variables: temperature (°C), salinity (ppt) and dissolved oxygen (mgL) and visibility (m).

Save the plot in tiff format

ggsave("Latitudevsabiotic.jpeg", units="mm", width=75, height=125, dpi=300)

Figure 3

Visibility_plot + Kelp_plot + Bare_plot + Inv_plot + PredFish_plot + Size_plot + Richness_plot + PCOA_plot+ Herb_plot + plot_layout(ncol=3)

Figure 3: Latitudinal variation of biotic and habitat variables

ggsave("FigurelatvsBio.jpeg", units="mm", width=150, height=100, dpi=300)

Phloro_plot

ggsave("phloro1.jpeg", units="mm", width=90, height=60, dpi=300)

Part 4: Latitudinal patterns in consumption

Latitudinal patterns in consumption (predation and herbivory) were explored and analysed

4a): Predation vs latitude

First change latitude to positive numbers for plotting

Predation_env$Latitude <- Predation_env$Latitude*-1

C_24 <- filter(Predation_env, Bait2 == "Crabs" & Assay_time2 == 24)
C_2 <- filter(Predation_env, Bait2 == "Crabs" & Assay_time2 == 2)
S_24 <- filter(Predation_env, Bait2 == "Squidpops" & Assay_time2 == 24)
S_2 <- filter(Predation_env, Bait2 == "Squidpops" & Assay_time2 == 2)

1) Crabs 2h

a. Analyses

# First, examine a simple scatter plot (always a good idea!)
   plot(C_2$Latitude, C_2$Consumed_individuals, xlab = "Latitude", ylab = "crabs consumed", las=1)
   mygrey <- grey(level = 0.65, alpha = 0.4)
   points(C_2$Latitude, C_2$Consumed_individuals, pch = 21, col = mygrey, bg = mygrey)

C2_mod1 <- mgcv::gam(Consumed_individuals ~ s(Latitude, k=11), gamma=1.4, data = C_2, family =poisson())
C2_mod2 <- mgcv::gam(Consumed_individuals ~ s(Latitude, k=11), gamma=1.4, data = C_2, family =nb())

AICc(C2_mod2, C2_mod1)

##                df     AICc
## C2_mod2  5.106187 446.0867
## C2_mod1 10.781434 428.1046

summary(C2_mod1)

## 
## Family: poisson 
## Link function: log 
## 
## Formula:
## Consumed_individuals ~ s(Latitude, k = 11)
## 
## Parametric coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   0.3941     0.0873   4.514 6.37e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##               edf Ref.df Chi.sq p-value    
## s(Latitude) 9.781  9.985  82.34  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.404   Deviance explained = 36.5%
## UBRE = 0.8541  Scale est. = 1         n = 120

par(mfrow = c(2, 2))
gam.check(C2_mod1)

## 
## Method: UBRE   Optimizer: outer newton
## full convergence after 6 iterations.
## Gradient range [4.749748e-07,4.749748e-07]
## (score 0.8540992 & scale 1).
## Hessian positive definite, eigenvalue range [0.003674509,0.003674509].
## Model rank =  11 / 11 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                k'   edf k-index p-value
## s(Latitude) 10.00  9.78    0.91    0.26

2) Crabs 24h

a. Analyses

# First, examine a simple scatter plot (always a good idea!)
   plot(C_24$Latitude, C_24$Consumed_individuals, xlab = "Latitude", ylab = "crabs consumed 24 h", las=1)
   mygrey <- grey(level = 0.65, alpha = 0.4)
   points(C_24$Latitude, C_24$Consumed_individuals, pch = 21, col = mygrey, bg = mygrey)

C24_mod1 <- mgcv::gam(Consumed_individuals ~ s(Latitude, k=11), gamma=1.4, data = C_24, family =poisson())
C24_mod2 <- mgcv::gam(Consumed_individuals ~ s(Latitude, k=11), gamma=1.4, data = C_24, family =nb())

AICc(C24_mod1, C24_mod2)

##                df     AICc
## C24_mod1 2.098819 433.6553
## C24_mod2 3.302821 436.0297

summary(C24_mod1)

## 
## Family: poisson 
## Link function: log 
## 
## Formula:
## Consumed_individuals ~ s(Latitude, k = 11)
## 
## Parametric coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.29862    0.04852   26.76   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##               edf Ref.df Chi.sq  p-value    
## s(Latitude) 1.099   1.19  29.66 9.54e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.385   Deviance explained = 34.2%
## UBRE = -0.46791  Scale est. = 1         n = 120

par(mfrow = c(2, 2))
gam.check(C24_mod1)

## 
## Method: UBRE   Optimizer: outer newton
## full convergence after 5 iterations.
## Gradient range [-3.566326e-08,-3.566326e-08]
## (score -0.467908 & scale 1).
## Hessian positive definite, eigenvalue range [0.0002048148,0.0002048148].
## Model rank =  11 / 11 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##               k'  edf k-index p-value  
## s(Latitude) 10.0  1.1    0.84    0.03 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Calculate predictions for Crab consumption after 2 and 24 h

#create a new dataframe: names should match those of the original dataframe
newdfr = expand.grid(Latitude = rep(seq(from = min(C_2$Latitude), to = max(C_2$Latitude), length.out = 100), 5))
 
#Crabs 2h

dat1 <- predict(C2_mod1, newdata = newdfr, type = "link", se.fit = TRUE)
critval <- 1.96 ## approx 95% CI
upr <- dat1$fit + (critval * dat1$se.fit)
lwr <- dat1$fit - (critval * dat1$se.fit)
fit <- dat1$fit
fit2 <- C2_mod1$family$linkinv(fit)
upr2 <- C2_mod1$family$linkinv(upr)
lwr2 <- C2_mod1$family$linkinv(lwr)

newdfr$lwr <- lwr2 
newdfr$upr <- upr2 
newdfr$fit <- fit2

#For crabs 24 h 

newdfr24 = expand.grid(Latitude = rep(seq(from = min(C_24$Latitude), to = max(C_24$Latitude), length.out = 100), 5))
 
dat2 <- predict(C24_mod1, newdata = newdfr24, type = "link", se.fit = TRUE)
upr1 <- dat2$fit + (critval * dat2$se.fit)
lwr1 <- dat2$fit - (critval * dat2$se.fit)
fit1 <- dat2$fit
fit3 <- C24_mod1$family$linkinv(fit1)
upr3 <- C24_mod1$family$linkinv(upr1)
lwr3 <- C24_mod1$family$linkinv(lwr1)

newdfr24$lwr <- lwr3 
newdfr24$upr <- upr3 
newdfr24$fit <- fit3

Add new column to mix datasets

newdfr$Assay_time <- "2 h"
newdfr24$Assay_time <- "24 h"
C_2$Assay_time <- "2 h"
C_24$Assay_time <- "24 h"

#Bind datasets

total <- rbind(newdfr, newdfr24)

total2 <- rbind (C_2, C_24)

3) Crabs 2 and 24 h prediction plot

ggplot()+
  geom_ribbon(data = total, aes(x = Latitude, y = fit, ymin = lwr, ymax = upr, group = Assay_time, fill = Assay_time ), alpha = .25)  +
  geom_line(data = total, aes(x = Latitude, y = fit, color = Assay_time, group = Assay_time), size = 1) +
  geom_jitter(data = total2, aes(x= Latitude, y=Consumed_individuals, color =Assay_time), width=0.15, height = 0.2, size = 1, alpha = 0.5)+
  labs(x = "Latitude", y = "Crabs consumed")  + Theme1 +
  scale_color_manual("Assay time",values=c("darkslateblue", "orange1"))+
   scale_fill_manual("Assay time", values=c("darkslateblue", "orange1")) + 
  theme(legend.position = c(0.88, 0.83),legend.title = element_text(size = 8), 
               legend.text = element_text(size = 8)) +
  scale_x_continuous(breaks = seq(15, 40, 5)) +
  scale_y_continuous(limits = c(0,6.5), breaks = seq(0, 5, 1))

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2
## 3.5.0.
## ℹ Please use the `legend.position.inside` argument of `theme()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

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

Save plot in tiff format

ggsave("Crabs.jpeg", units="mm", width=120, height=75, dpi=300)

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

4) Squidpops 2h

a. Analyses

S2_mod1 <- mgcv::gam(Consumed_individuals ~ s(Latitude, k=10), gamma=1.4, data = S_2, family =poisson())
S2_mod2 <- mgcv::gam(Consumed_individuals ~ s(Latitude, k=10), gamma=1.4, S_2, family =nb())

AICc(S2_mod2, S2_mod1)

##               df     AICc
## S2_mod2 3.001418 490.8237
## S2_mod1 9.825399 462.7014

summary(S2_mod1)

## 
## Family: poisson 
## Link function: log 
## 
## Formula:
## Consumed_individuals ~ s(Latitude, k = 10)
## 
## Parametric coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.72206    0.06967   10.36   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##               edf Ref.df Chi.sq p-value    
## s(Latitude) 8.825  8.989  50.88  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.209   Deviance explained = 27.2%
## UBRE = 0.77145  Scale est. = 1         n = 120

par(mfrow = c(2, 2))
gam.check(S2_mod1)

## 
## Method: UBRE   Optimizer: outer newton
## full convergence after 6 iterations.
## Gradient range [1.464168e-06,1.464168e-06]
## (score 0.7714456 & scale 1).
## Hessian positive definite, eigenvalue range [0.003165155,0.003165155].
## Model rank =  10 / 10 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##               k'  edf k-index p-value    
## s(Latitude) 9.00 8.83    0.56  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

5) Squidpops 24h

a. Analyses

S24_mod1 <- mgcv::gam(Consumed_individuals ~ s(Latitude, k=10), gamma=1.4, data = S_24, family =poisson())
S24_mod2 <- mgcv::gam(Consumed_individuals ~ s(Latitude, k=10), gamma=1.4, S_24, family =nb())

AICc(S24_mod2, S24_mod1)

##                df     AICc
## S24_mod2 3.000007 448.2412
## S24_mod1 2.000015 446.1367

summary(S24_mod1)

## 
## Family: poisson 
## Link function: log 
## 
## Formula:
## Consumed_individuals ~ s(Latitude, k = 10)
## 
## Parametric coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.43953    0.04461   32.27   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##             edf Ref.df Chi.sq p-value   
## s(Latitude)   1      1   7.33 0.00678 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.154   Deviance explained = 12.1%
## UBRE = -0.501  Scale est. = 1         n = 120

par(mfrow = c(2, 2))
gam.check(S24_mod1)

## 
## Method: UBRE   Optimizer: outer newton
## full convergence after 10 iterations.
## Gradient range [-2.01421e-07,-2.01421e-07]
## (score -0.5009978 & scale 1).
## Hessian positive definite, eigenvalue range [2.014182e-07,2.014182e-07].
## Model rank =  10 / 10 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##             k' edf k-index p-value    
## s(Latitude)  9   1    0.75  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#For squidpops 24 h 

newdata_S24 <- data.frame(Latitude = rep(seq(from = min(S_24$Latitude), to = max(S_24$Latitude), length.out = 100), 5))

#add preditions
dat3 <- predict(S24_mod1, newdata = newdata_S24, type = "link", se.fit = TRUE)
upr1 <- dat3$fit + (critval * dat3$se.fit)
lwr1<- dat3$fit - (critval * dat3$se.fit)
fit1 <- dat3$fit
fit2 <- S24_mod1$family$linkinv(fit1)
upr2 <- S24_mod1$family$linkinv(upr1)
lwr2 <- S24_mod1$family$linkinv(lwr1)

newdata_S24$lwr <- lwr2 
newdata_S24$upr <- upr2 
newdata_S24$fit <- fit2

#For squidpops 24 h 

newdata_S2 <- data.frame(Latitude = rep(seq(from = min(S_2$Latitude), to = max(S_2$Latitude), length.out = 100), 5))

#add preditions
dat4 <- predict(S2_mod1, newdata = newdata_S24, type = "link", se.fit = TRUE)
upr4 <- dat4$fit + (critval * dat4$se.fit)
lwr4<- dat4$fit - (critval * dat4$se.fit)
fit4 <- dat4$fit
fit5 <- S2_mod1$family$linkinv(fit4)
upr5 <- S2_mod1$family$linkinv(upr4)
lwr5 <- S2_mod1$family$linkinv(lwr4)

newdata_S2$lwr <- lwr5 
newdata_S2$upr <- upr5 
newdata_S2$fit <- fit5

Add new column to mix datasets

S_2$Assay_time <- "2 h"
S_24$Assay_time <- "24 h"
newdata_S24$Assay_time <- "24 h"
newdata_S2$Assay_time <- "2 h"
##Keep only colums for plotting
total <- rbind (newdata_S2, newdata_S24)
total2 <- rbind (S_2, S_24)

3) Squidpops 2 and 24 h prediction plot

ggplot() +
  geom_ribbon(data = total, aes(x = Latitude, y = fit, ymin = lwr, ymax = upr, group = Assay_time, fill = Assay_time ), alpha = .25) +
  geom_line(data = total, aes(x = Latitude, y = fit, color = Assay_time, group = Assay_time ), size = 1) +
  geom_jitter(data = total2, aes(x= Latitude, y=Consumed_individuals, color =Assay_time), width=0.15, height = 0.2, size = 1, alpha = .5)+
  labs(x = "Latitude", y = "Squidpops consumed")  + Theme1 +
  scale_color_manual("Assay time",values=c("darkslateblue", "orange1"))+
   scale_fill_manual("Assay time", values=c("darkslateblue", "orange1")) + 
  theme(legend.position = "none") +
  scale_x_continuous(breaks = seq(15, 40, 5)) +
  scale_y_continuous(limits = c(0,6.5), breaks = seq(0, 5, 1))

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

Save plot in tiff format

ggsave("Squidpops.jpeg", units="mm", width=120, height=70, dpi=300)

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

4b): Herbivory vs latitude

Herbivory_env$Latitude <- Herbivory_env$Latitude*-1

a. Analyses

Herb_mod1 <- mgcv::gam(Consumed_proportion ~ s(Latitude), data = Herbivory_env, family =binomial())

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

Herb_mod2 <- mgcv::gam(Consumed_proportion ~ s(Latitude), data = Herbivory_env, family =quasibinomial())

AICc(Herb_mod2, Herb_mod1)

##                  df     AICc
## Herb_mod2 10.863328       NA
## Herb_mod1  9.759176 413.4238

summary(Herb_mod1)

## 
## Family: binomial 
## Link function: logit 
## 
## Formula:
## Consumed_proportion ~ s(Latitude)
## 
## Parametric coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -3.1503     0.5485  -5.743 9.29e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##               edf Ref.df Chi.sq p-value    
## s(Latitude) 8.759  8.977  81.65  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.451   Deviance explained = 46.4%
## UBRE = -0.52027  Scale est. = 1         n = 590

par(mfrow = c(2, 2))
gam.check(Herb_mod1)

## 
## Method: UBRE   Optimizer: outer newton
## full convergence after 8 iterations.
## Gradient range [5.631959e-09,5.631959e-09]
## (score -0.5202735 & scale 1).
## Hessian positive definite, eigenvalue range [0.000357682,0.000357682].
## Model rank =  10 / 10 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##               k'  edf k-index p-value    
## s(Latitude) 9.00 8.76    0.46  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#For squidpops 24 h 

newdata_H <- data.frame(Latitude = rep(seq(from = min(Herbivory_env$Latitude), to = max(Herbivory_env$Latitude), length.out = 100), 5))

#add preditions
dat3 <- predict(Herb_mod1, newdata = newdata_H, type = "link", se.fit = TRUE)
upr1 <- dat3$fit + (critval * dat3$se.fit)
lwr1<- dat3$fit - (critval * dat3$se.fit)
fit1 <- dat3$fit
fit2 <- Herb_mod1$family$linkinv(fit1)
upr2 <- Herb_mod1$family$linkinv(upr1)
lwr2 <- Herb_mod1$family$linkinv(lwr1)

newdata_H$lwr <- lwr2 
newdata_H$upr <- upr2 
newdata_H$fit <- fit2

b. Plot predicted relationship

ggplot() +
  geom_ribbon(data = newdata_H, aes(x = Latitude, y = fit, ymin = lwr, ymax = upr), alpha = .25, fill = "orange1") +
  geom_line(data = newdata_H, aes(x = Latitude, y = fit ), size = 1, color = "orange1") +
  geom_jitter(data = Herbivory_env, aes(x= Latitude, y=Consumed_proportion), color = "orange1", width=0.2, height = 0.02, size = 1, alpha = .5)+
  labs(x = "Latitude (°S)", y = "Blade consumed proportion")  + Theme1  + 
  theme(legend.position = "none") +
  scale_x_continuous(breaks = seq(15, 40, 5)) +
  scale_y_continuous(limits = c(-0.1,1.2), breaks = seq(0, 1, 0.2))

Save plot to tiff format

ggsave("Seaweed24.jpeg", units="mm", width=120, height=75, dpi=300)

Part 5: Predation intensity: Effect of site-specific variables

Here, the individual effects of each predictor will be tested.

First, create datasets to be used in the analysis

Predation_analysis <- Predation_6 %>%
  dplyr::select(-Rope, -Bait, -Camera, -Assay_time, -Initial_Individuals, -Remaining_individuals, -Survival_proportion,-Consumed_proportion)

#Separate by Bait and Assay time 

Crabs_24 <- filter(Predation_analysis, Bait2 == "Crabs" & Assay_time2 == 24)
Crabs_2 <- filter(Predation_analysis, Bait2 == "Crabs" & Assay_time2 == 2)
Squidpops_24 <- filter(Predation_analysis, Bait2 == "Squidpops" & Assay_time2 == 24)
Squidpops_2 <- filter(Predation_analysis, Bait2 == "Squidpops" & Assay_time2 == 2)

5a): Standardisation and test autocorrelation

1) Standardise

#Crabs 24h
Crabs24_num <- Crabs_24 %>%
  dplyr::select(-Site, -Bait2, -Assay_time2, -Consumed_individuals, -Latitude, -Day_length) 
Crabs24_stand <- decostand(Crabs24_num, method="range", na.rm = TRUE)
#Test autocorrelation 
corrplot(cor(Crabs24_stand), method = "number")

#Add site and consumed individuals for analyses
Crabs24_stand$Site <- Crabs_24$Site
Crabs24_stand$Consumed_individuals<- Crabs_24$Consumed_individuals
Crabs24_stand$Latitude<- Crabs_24$Latitude

No multicolinearity found

Standardise all final variables to be in a range 0-1

#Crabs 2h 
Crabs2_num <- Crabs_2 %>%
  dplyr::select(-Site, -Bait2, -Assay_time2, -Consumed_individuals, -Latitude) 
Crabs2_stand <- decostand(Crabs2_num, method="range", na.rm = TRUE)
Crabs2_stand$Site <- Crabs_2$Site
Crabs2_stand$Consumed_individuals<- Crabs_2$Consumed_individuals
Crabs2_stand$Latitude<- Crabs_2$Latitude
#Squidpops 24h
Squidpops24_num <- Squidpops_24 %>%
  dplyr::select(-Site, -Bait2, -Assay_time2, -Consumed_individuals, -Latitude) 
Squidpops24_stand <- decostand(Squidpops24_num, method="range", na.rm = TRUE)
Squidpops24_stand$Site <- Squidpops_24$Site
Squidpops24_stand$Consumed_individuals<- Squidpops_24$Consumed_individuals
Squidpops24_stand$Latitude<- Squidpops_24$Latitude
#Squidpops 2h
Squidpops2_num <- Squidpops_2 %>%
  dplyr::select(-Site, -Bait2, -Assay_time2, -Consumed_individuals, -Latitude) 
Squidpops2_stand <- decostand(Squidpops2_num, method="range", na.rm = TRUE)
Squidpops2_stand$Site <- Squidpops_2$Site
Squidpops2_stand$Consumed_individuals<- Squidpops_2$Consumed_individuals
Squidpops2_stand$Latitude<- Squidpops_2$Latitude

5b) Model selection: Crabs 2h

Temperature

my_glmT<- glm(Consumed_individuals~  Temperature, data=Crabs2_stand, family = poisson)
check_overdispersion(my_glmT)

## # Overdispersion test
## 
##        dispersion ratio =   1.734
##   Pearson's Chi-Squared = 204.668
##                 p-value = < 0.001

## Overdispersion detected.

my_glmT<- glm.nb(Consumed_individuals~  Temperature, data=Crabs2_stand)
check_overdispersion(my_glmT)

## # Overdispersion test
## 
##  dispersion ratio = 0.666
##           p-value =  0.08

## No overdispersion detected.

summary(my_glmT)

## 
## Call:
## glm.nb(formula = Consumed_individuals ~ Temperature, data = Crabs2_stand, 
##     init.theta = 1.804741369, link = log)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.7265     0.2970  -2.446   0.0144 *  
## Temperature   1.9321     0.3861   5.003 5.63e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(1.8047) family taken to be 1)
## 
##     Null deviance: 173.79  on 119  degrees of freedom
## Residual deviance: 143.22  on 118  degrees of freedom
## AIC: 435.29
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  1.805 
##           Std. Err.:  0.633 
## 
##  2 x log-likelihood:  -429.289

AICc(my_glmT)

## [1] 435.4963

Salinity

my_glmSC2<- glm(Consumed_individuals~  Salinity_ppt, data=Crabs2_stand, family = poisson)
check_overdispersion(my_glmSC2)

## # Overdispersion test
## 
##        dispersion ratio =   2.003
##   Pearson's Chi-Squared = 236.411
##                 p-value = < 0.001

## Overdispersion detected.

my_glmSC2<- glm.nb(Consumed_individuals~  Salinity_ppt, data=Crabs2_stand)
check_overdispersion(my_glmSC2)

## # Overdispersion test
## 
##  dispersion ratio = 0.702
##           p-value = 0.184

## No overdispersion detected.

summary(my_glmSC2)

## 
## Call:
## glm.nb(formula = Consumed_individuals ~ Salinity_ppt, data = Crabs2_stand, 
##     init.theta = 1.09232541, link = log)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)   
## (Intercept)   -0.4038     0.3777  -1.069  0.28507   
## Salinity_ppt   1.5256     0.5117   2.982  0.00287 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(1.0923) family taken to be 1)
## 
##     Null deviance: 141.19  on 119  degrees of freedom
## Residual deviance: 133.24  on 118  degrees of freedom
## AIC: 452.5
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  1.092 
##           Std. Err.:  0.282 
## 
##  2 x log-likelihood:  -446.499

AICc(my_glmSC2)

## [1] 452.7059

Dissolved oxygen

my_glmOx<- glm.nb(Consumed_individuals~  D_oxygen_mgL, data=Crabs2_stand)
summary(my_glmOx)

## 
## Call:
## glm.nb(formula = Consumed_individuals ~ D_oxygen_mgL, data = Crabs2_stand, 
##     init.theta = 1.455072917, link = log)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -0.0611     0.1946  -0.314    0.754    
## D_oxygen_mgL   1.4488     0.3371   4.298 1.72e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(1.4551) family taken to be 1)
## 
##     Null deviance: 159.78  on 119  degrees of freedom
## Residual deviance: 139.58  on 118  degrees of freedom
## AIC: 442.61
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  1.455 
##           Std. Err.:  0.452 
## 
##  2 x log-likelihood:  -436.609

AICc(my_glmOx)

## [1] 442.8161

Visibility

my_glm<- glm.nb(Consumed_individuals~ V_Visibility, data=Crabs2_stand)
summary(my_glm)

## 
## Call:
## glm.nb(formula = Consumed_individuals ~ V_Visibility, data = Crabs2_stand, 
##     init.theta = 0.9594356772, link = log)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)   
## (Intercept)    0.5989     0.2140   2.799  0.00512 **
## V_Visibility   0.1710     0.4245   0.403  0.68703   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(0.9594) family taken to be 1)
## 
##     Null deviance: 132.93  on 119  degrees of freedom
## Residual deviance: 132.79  on 118  degrees of freedom
## AIC: 460.03
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  0.959 
##           Std. Err.:  0.235 
## 
##  2 x log-likelihood:  -454.031

AICc(my_glm)

## [1] 460.2381

Community composition

my_glm<- glm.nb(Consumed_individuals~  PCo1, data=Crabs2_stand)
summary(my_glm)

## 
## Call:
## glm.nb(formula = Consumed_individuals ~ PCo1, data = Crabs2_stand, 
##     init.theta = 0.9572416306, link = log)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.68732    0.16118   4.264 2.01e-05 ***
## PCo1        -0.04288    0.32097  -0.134    0.894    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(0.9572) family taken to be 1)
## 
##     Null deviance: 132.79  on 119  degrees of freedom
## Residual deviance: 132.77  on 118  degrees of freedom
## AIC: 460.15
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  0.957 
##           Std. Err.:  0.235 
## 
##  2 x log-likelihood:  -454.154

AICc(my_glm)

## [1] 460.3613

Kelp Cover

my_glm<- glm.nb(Consumed_individuals~ Kelp, data=Crabs2_stand)
summary(my_glm)

## 
## Call:
## glm.nb(formula = Consumed_individuals ~ Kelp, data = Crabs2_stand, 
##     init.theta = 1.175060865, link = log)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)   0.3118     0.1547   2.015  0.04389 * 
## Kelp          1.1196     0.3639   3.077  0.00209 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(1.1751) family taken to be 1)
## 
##     Null deviance: 145.89  on 119  degrees of freedom
## Residual deviance: 135.96  on 118  degrees of freedom
## AIC: 450.9
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  1.175 
##           Std. Err.:  0.321 
## 
##  2 x log-likelihood:  -444.902

AICc(my_glm)

## [1] 451.1094

Bare rock cover

my_glmBR<- glm.nb(Consumed_individuals~ Bare_rock, data=Crabs2_stand)
summary(my_glmBR)

## 
## Call:
## glm.nb(formula = Consumed_individuals ~ Bare_rock, data = Crabs2_stand, 
##     init.theta = 1.228290388, link = log)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   1.1046     0.1536   7.190 6.48e-13 ***
## Bare_rock    -1.3545     0.3449  -3.927 8.58e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(1.2283) family taken to be 1)
## 
##     Null deviance: 148.76  on 119  degrees of freedom
## Residual deviance: 134.51  on 118  degrees of freedom
## AIC: 446.9
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  1.228 
##           Std. Err.:  0.335 
## 
##  2 x log-likelihood:  -440.900

AICc(my_glmBR)

## [1] 447.1067

Species richness

my_glm<- glm.nb(Consumed_individuals~ Mean_SR, data=Crabs2_stand)
summary(my_glm)

## 
## Call:
## glm.nb(formula = Consumed_individuals ~ Mean_SR, data = Crabs2_stand, 
##     init.theta = 1.058830107, link = log)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)   0.2469     0.2247   1.099   0.2719  
## Mean_SR       0.8301     0.3910   2.123   0.0337 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(1.0588) family taken to be 1)
## 
##     Null deviance: 139.20  on 119  degrees of freedom
## Residual deviance: 134.34  on 118  degrees of freedom
## AIC: 455.48
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  1.059 
##           Std. Err.:  0.274 
## 
##  2 x log-likelihood:  -449.479

AICc(my_glm)

## [1] 455.6857

Invertebrate predators

my_glm<- glm.nb(Consumed_individuals~ Invertebrate_predators, data=Crabs2_stand)
summary(my_glm)

## 
## Call:
## glm.nb(formula = Consumed_individuals ~ Invertebrate_predators, 
##     data = Crabs2_stand, init.theta = 1.184638797, link = log)
## 
## Coefficients:
##                        Estimate Std. Error z value Pr(>|z|)   
## (Intercept)              0.3891     0.1406   2.768  0.00564 **
## Invertebrate_predators   0.8476     0.2939   2.884  0.00393 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(1.1846) family taken to be 1)
## 
##     Null deviance: 146.42  on 119  degrees of freedom
## Residual deviance: 136.66  on 118  degrees of freedom
## AIC: 451.13
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  1.185 
##           Std. Err.:  0.328 
## 
##  2 x log-likelihood:  -445.133

AICc(my_glm)

## [1] 451.3399

Predatory fish abundance

my_glm <- glm.nb(Consumed_individuals~ Predatory_fish, data=Crabs2_stand)
summary(my_glm)

## 
## Call:
## glm.nb(formula = Consumed_individuals ~ Predatory_fish, data = Crabs2_stand, 
##     init.theta = 0.975580135, link = log)
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)   
## (Intercept)      0.5454     0.1713   3.183  0.00146 **
## Predatory_fish   0.3192     0.3293   0.969  0.33239   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(0.9756) family taken to be 1)
## 
##     Null deviance: 133.99  on 119  degrees of freedom
## Residual deviance: 133.04  on 118  degrees of freedom
## AIC: 459.23
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  0.976 
##           Std. Err.:  0.241 
## 
##  2 x log-likelihood:  -453.226

AICc(my_glm)

## [1] 459.4329

Fish size

my_glm<- glm.nb(Consumed_individuals~ Fish_size, data=Crabs2_stand)
summary(my_glm)

## 
## Call:
## glm.nb(formula = Consumed_individuals ~ Fish_size, data = Crabs2_stand, 
##     init.theta = 0.967006894, link = log)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)   0.3846     0.2900   1.326    0.185
## Fish_size     0.4582     0.4246   1.079    0.281
## 
## (Dispersion parameter for Negative Binomial(0.967) family taken to be 1)
## 
##     Null deviance: 133.43  on 119  degrees of freedom
## Residual deviance: 132.69  on 118  degrees of freedom
## AIC: 459.43
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  0.967 
##           Std. Err.:  0.237 
## 
##  2 x log-likelihood:  -453.429

AICc(my_glm)

## [1] 459.6363

Test interaction between temperature and dissolved oxygen

my_glm<- glm.nb(Consumed_individuals~ Temperature*D_oxygen_mgL, data=Crabs2_stand)
summary(my_glm)

## 
## Call:
## glm.nb(formula = Consumed_individuals ~ Temperature * D_oxygen_mgL, 
##     data = Crabs2_stand, init.theta = 2.355100978, link = log)
## 
## Coefficients:
##                          Estimate Std. Error z value Pr(>|z|)  
## (Intercept)              -0.55352    0.42946  -1.289   0.1974  
## Temperature               1.15886    0.55926   2.072   0.0383 *
## D_oxygen_mgL             -0.09859    1.13581  -0.087   0.9308  
## Temperature:D_oxygen_mgL  1.08330    1.25059   0.866   0.3864  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(2.3551) family taken to be 1)
## 
##     Null deviance: 190.78  on 119  degrees of freedom
## Residual deviance: 148.41  on 116  degrees of freedom
## AIC: 432.35
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  2.355 
##           Std. Err.:  0.998 
## 
##  2 x log-likelihood:  -422.350

AICc(my_glm)

## [1] 432.8764

5c) Model selection: Crabs 24h

Temperature

my_glmT24<- glm(Consumed_individuals~  Temperature, data=Crabs24_stand, family = poisson)
check_overdispersion(my_glmT24)

## # Overdispersion test
## 
##        dispersion ratio =  0.431
##   Pearson's Chi-Squared = 50.817
##                 p-value =      1

## No overdispersion detected.

summary(my_glmT24)

## 
## Call:
## glm(formula = Consumed_individuals ~ Temperature, family = poisson, 
##     data = Crabs24_stand)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   0.7381     0.1345   5.487 4.09e-08 ***
## Temperature   0.8747     0.1767   4.950 7.41e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 88.137  on 119  degrees of freedom
## Residual deviance: 61.359  on 118  degrees of freedom
## AIC: 436.73
## 
## Number of Fisher Scoring iterations: 4

AICc(my_glmT24)

## [1] 436.8334

Salinity

my_glmS24<- glm(Consumed_individuals~  Salinity_ppt, data=Crabs24_stand, family = poisson)
summary(my_glmS24)

## 
## Call:
## glm(formula = Consumed_individuals ~ Salinity_ppt, family = poisson, 
##     data = Crabs24_stand)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    0.6381     0.1695   3.765 0.000167 ***
## Salinity_ppt   0.9948     0.2255   4.412 1.02e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 88.137  on 119  degrees of freedom
## Residual deviance: 66.419  on 118  degrees of freedom
## AIC: 441.79
## 
## Number of Fisher Scoring iterations: 4

AICc(my_glmS24)

## [1] 441.8936

Dissolved oxygen

my_glm<- glm(Consumed_individuals~  D_oxygen_mgL, data=Crabs24_stand, family = poisson)
summary(my_glm)

## 
## Call:
## glm(formula = Consumed_individuals ~ D_oxygen_mgL, family = poisson, 
##     data = Crabs24_stand)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   1.22985    0.08573  14.345   <2e-16 ***
## D_oxygen_mgL  0.23004    0.15682   1.467    0.142    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 88.137  on 119  degrees of freedom
## Residual deviance: 86.014  on 118  degrees of freedom
## AIC: 461.39
## 
## Number of Fisher Scoring iterations: 4

AICc(my_glm)

## [1] 461.4883

Visibility

my_glm<- glm(Consumed_individuals~ V_Visibility, data=Crabs24_stand, family = poisson)
summary(my_glm)

## 
## Call:
## glm(formula = Consumed_individuals ~ V_Visibility, family = poisson, 
##     data = Crabs24_stand)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   1.43414    0.08672  16.537   <2e-16 ***
## V_Visibility -0.24462    0.17978  -1.361    0.174    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 88.137  on 119  degrees of freedom
## Residual deviance: 86.257  on 118  degrees of freedom
## AIC: 461.63
## 
## Number of Fisher Scoring iterations: 4

AICc(my_glm)

## [1] 461.7317

Consumer community

my_glm<- glm(Consumed_individuals~  PCo1, data=Crabs24_stand, family = poisson)
summary(my_glm)

## 
## Call:
## glm(formula = Consumed_individuals ~ PCo1, family = poisson, 
##     data = Crabs24_stand)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.45277    0.06421  22.625   <2e-16 ***
## PCo1        -0.35762    0.13953  -2.563   0.0104 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 88.137  on 119  degrees of freedom
## Residual deviance: 81.302  on 118  degrees of freedom
## AIC: 456.67
## 
## Number of Fisher Scoring iterations: 4

AICc(my_glm)

## [1] 456.7768

Kelp cover

my_glm<- glm(Consumed_individuals~ Kelp, data=Crabs24_stand, family = poisson)
summary(my_glm)

## 
## Call:
## glm(formula = Consumed_individuals ~ Kelp, family = poisson, 
##     data = Crabs24_stand)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.17237    0.06767  17.326  < 2e-16 ***
## Kelp         0.54234    0.15254   3.555 0.000377 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 88.137  on 119  degrees of freedom
## Residual deviance: 76.268  on 118  degrees of freedom
## AIC: 451.64
## 
## Number of Fisher Scoring iterations: 4

AICc(my_glm)

## [1] 451.7428

Bare rock cover

my_glmBR24<- glm(Consumed_individuals~ Bare_rock, data=Crabs24_stand, family = poisson)
summary(my_glmBR24)

## 
## Call:
## glm(formula = Consumed_individuals ~ Bare_rock, family = poisson, 
##     data = Crabs24_stand)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.52718    0.06716  22.740  < 2e-16 ***
## Bare_rock   -0.55311    0.14903  -3.711 0.000206 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 88.137  on 119  degrees of freedom
## Residual deviance: 73.666  on 118  degrees of freedom
## AIC: 449.04
## 
## Number of Fisher Scoring iterations: 4

AICc(my_glmBR24)

## [1] 449.1402

Species richness

my_glm<- glm(Consumed_individuals~ Mean_SR, data=Crabs24_stand, family = poisson)
summary(my_glm)

## 
## Call:
## glm(formula = Consumed_individuals ~ Mean_SR, family = poisson, 
##     data = Crabs24_stand)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.15922    0.09538  12.154   <2e-16 ***
## Mean_SR      0.35367    0.16445   2.151   0.0315 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 88.137  on 119  degrees of freedom
## Residual deviance: 83.537  on 118  degrees of freedom
## AIC: 458.91
## 
## Number of Fisher Scoring iterations: 4

AICc(my_glm)

## [1] 459.0114

Abundance of invertebrate predators

my_glm<- glm(Consumed_individuals~ Invertebrate_predators, data=Crabs24_stand, family = poisson)
summary(my_glm)

## 
## Call:
## glm(formula = Consumed_individuals ~ Invertebrate_predators, 
##     family = poisson, data = Crabs24_stand)
## 
## Coefficients:
##                        Estimate Std. Error z value Pr(>|z|)    
## (Intercept)              1.2383     0.0614  20.167   <2e-16 ***
## Invertebrate_predators   0.3246     0.1270   2.556   0.0106 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 88.137  on 119  degrees of freedom
## Residual deviance: 81.890  on 118  degrees of freedom
## AIC: 457.26
## 
## Number of Fisher Scoring iterations: 4

AICc(my_glm)

## [1] 457.3644

Abundance of fish predators

my_glm<- glm(Consumed_individuals~ Predatory_fish, data=Crabs24_stand, family = poisson)
summary(my_glm)

## 
## Call:
## glm(formula = Consumed_individuals ~ Predatory_fish, family = poisson, 
##     data = Crabs24_stand)
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     1.28137    0.07112  18.018   <2e-16 ***
## Predatory_fish  0.13325    0.13578   0.981    0.326    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 88.137  on 119  degrees of freedom
## Residual deviance: 87.184  on 118  degrees of freedom
## AIC: 462.56
## 
## Number of Fisher Scoring iterations: 4

AICc(my_glm)

## [1] 462.6586

Fish legth (cm)

my_glm24_FL<- glm(Consumed_individuals~ Fish_size, data=Crabs24_stand, family = poisson)
summary(my_glm24_FL)

## 
## Call:
## glm(formula = Consumed_individuals ~ Fish_size, family = poisson, 
##     data = Crabs24_stand)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   0.9420     0.1297   7.263  3.8e-13 ***
## Fish_size     0.6118     0.1833   3.337 0.000847 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 88.137  on 119  degrees of freedom
## Residual deviance: 76.519  on 118  degrees of freedom
## AIC: 451.89
## 
## Number of Fisher Scoring iterations: 4

AICc(my_glm24_FL)

## [1] 451.9938

Test interaction between temperature and oxygen

my_glm<- glm(Consumed_individuals~ Temperature*D_oxygen_mgL , data=Crabs24_stand, family = poisson)
summary(my_glm)

## 
## Call:
## glm(formula = Consumed_individuals ~ Temperature * D_oxygen_mgL, 
##     family = poisson, data = Crabs24_stand)
## 
## Coefficients:
##                          Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                0.9071     0.1883   4.817 1.46e-06 ***
## Temperature                0.7731     0.2518   3.070  0.00214 ** 
## D_oxygen_mgL              -0.7182     0.5435  -1.321  0.18641    
## Temperature:D_oxygen_mgL   0.6468     0.6046   1.070  0.28471    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 88.137  on 119  degrees of freedom
## Residual deviance: 59.281  on 116  degrees of freedom
## AIC: 438.65
## 
## Number of Fisher Scoring iterations: 4

AICc(my_glm)

## [1] 439.0009

5e): Model selection: Squidpops 2h

Temperature

my_glm<- glm.nb(Consumed_individuals~  Temperature, data=Squidpops2_stand)
summary(my_glm)

## 
## Call:
## glm.nb(formula = Consumed_individuals ~ Temperature, data = Squidpops2_stand, 
##     init.theta = 2.246867963, link = log)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)   0.5822     0.2111   2.758  0.00582 **
## Temperature   0.4259     0.2918   1.460  0.14435   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(2.2469) family taken to be 1)
## 
##     Null deviance: 150.71  on 119  degrees of freedom
## Residual deviance: 148.34  on 118  degrees of freedom
## AIC: 488.36
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  2.247 
##           Std. Err.:  0.719 
## 
##  2 x log-likelihood:  -482.357

AICc(my_glm)

## [1] 488.5636

Salinity

my_glm<- glm.nb(Consumed_individuals~  Salinity_ppt, data=Squidpops2_stand)
summary(my_glm)

## 
## Call:
## glm.nb(formula = Consumed_individuals ~ Salinity_ppt, data = Squidpops2_stand, 
##     init.theta = 2.124515491, link = log)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    0.9917     0.2464   4.025 5.69e-05 ***
## Salinity_ppt  -0.1904     0.3465  -0.550    0.583    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(2.1245) family taken to be 1)
## 
##     Null deviance: 147.65  on 119  degrees of freedom
## Residual deviance: 147.34  on 118  degrees of freedom
## AIC: 490.37
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  2.125 
##           Std. Err.:  0.655 
## 
##  2 x log-likelihood:  -484.369

AICc(my_glm)

## [1] 490.5761

Dissolved oxygen

my_glmOxS<- glm.nb(Consumed_individuals~  D_oxygen_mgL, data=Squidpops2_stand)
summary(my_glmOxS)

## 
## Call:
## glm.nb(formula = Consumed_individuals ~ D_oxygen_mgL, data = Squidpops2_stand, 
##     init.theta = 2.926891485, link = log)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    0.4251     0.1509   2.817 0.004844 ** 
## D_oxygen_mgL   0.9177     0.2651   3.462 0.000537 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(2.9269) family taken to be 1)
## 
##     Null deviance: 164.83  on 119  degrees of freedom
## Residual deviance: 152.57  on 118  degrees of freedom
## AIC: 479.42
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  2.93 
##           Std. Err.:  1.11 
## 
##  2 x log-likelihood:  -473.423

AICc(my_glmOxS)

## [1] 479.6298

Visibility

my_glm<- glm.nb(Consumed_individuals~ V_Visibility, data=Squidpops2_stand)
summary(my_glm)

## 
## Call:
## glm.nb(formula = Consumed_individuals ~ V_Visibility, data = Squidpops2_stand, 
##     init.theta = 2.191189762, link = log)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    0.6935     0.1623   4.272 1.94e-05 ***
## V_Visibility   0.3923     0.3161   1.241    0.215    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(2.1912) family taken to be 1)
## 
##     Null deviance: 149.34  on 119  degrees of freedom
## Residual deviance: 147.80  on 118  degrees of freedom
## AIC: 489.15
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  2.191 
##           Std. Err.:  0.688 
## 
##  2 x log-likelihood:  -483.153

AICc(my_glm)

## [1] 489.3599

Community composition

my_glmPC<- glm.nb(Consumed_individuals~  PCo1, data=Squidpops2_stand)
summary(my_glmPC)

## 
## Call:
## glm.nb(formula = Consumed_individuals ~ PCo1, data = Squidpops2_stand, 
##     init.theta = 2.570655523, link = log)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   0.5998     0.1221   4.913 8.98e-07 ***
## PCo1          0.6659     0.2237   2.977  0.00291 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(2.5707) family taken to be 1)
## 
##     Null deviance: 157.98  on 119  degrees of freedom
## Residual deviance: 149.71  on 118  degrees of freedom
## AIC: 482.8
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  2.571 
##           Std. Err.:  0.887 
## 
##  2 x log-likelihood:  -476.801

AICc(my_glmPC)

## [1] 483.0077

Kelp cover

my_glm<- glm.nb(Consumed_individuals~ Kelp, data=Squidpops2_stand)
summary(my_glm)

## 
## Call:
## glm.nb(formula = Consumed_individuals ~ Kelp, data = Squidpops2_stand, 
##     init.theta = 2.435059357, link = log)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   0.6668     0.1193   5.591 2.26e-08 ***
## Kelp          0.6583     0.2841   2.317   0.0205 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(2.4351) family taken to be 1)
## 
##     Null deviance: 155.07  on 119  degrees of freedom
## Residual deviance: 149.59  on 118  degrees of freedom
## AIC: 485.42
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  2.435 
##           Std. Err.:  0.820 
## 
##  2 x log-likelihood:  -479.418

AICc(my_glm)

## [1] 485.6245

Bare rock cover

my_glm<- glm.nb(Consumed_individuals~ Bare_rock, data=Squidpops2_stand)
summary(my_glm)

## 
## Call:
## glm.nb(formula = Consumed_individuals ~ Bare_rock, data = Squidpops2_stand, 
##     init.theta = 2.218274096, link = log)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   0.9933     0.1267   7.838 4.59e-15 ***
## Bare_rock    -0.3555     0.2587  -1.374    0.169    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(2.2183) family taken to be 1)
## 
##     Null deviance: 150.01  on 119  degrees of freedom
## Residual deviance: 148.03  on 118  degrees of freedom
## AIC: 488.73
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  2.218 
##           Std. Err.:  0.703 
## 
##  2 x log-likelihood:  -482.727

AICc(my_glm)

## [1] 488.9339

Species richness

my_glm<- glm.nb(Consumed_individuals~ Mean_SR, data=Squidpops2_stand)
summary(my_glm)

## 
## Call:
## glm.nb(formula = Consumed_individuals ~ Mean_SR, data = Squidpops2_stand, 
##     init.theta = 2.108078559, link = log)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.84656    0.16972   4.988 6.11e-07 ***
## Mean_SR      0.03856    0.30561   0.126      0.9    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(2.1081) family taken to be 1)
## 
##     Null deviance: 147.23  on 119  degrees of freedom
## Residual deviance: 147.21  on 118  degrees of freedom
## AIC: 490.67
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  2.108 
##           Std. Err.:  0.647 
## 
##  2 x log-likelihood:  -484.667

AICc(my_glm)

## [1] 490.8737

Abundance of invertebrate predators

my_glm<- glm.nb(Consumed_individuals~ Invertebrate_predators, data=Squidpops2_stand)
summary(my_glm)

## 
## Call:
## glm.nb(formula = Consumed_individuals ~ Invertebrate_predators, 
##     data = Squidpops2_stand, init.theta = 2.21273312, link = log)
## 
## Coefficients:
##                        Estimate Std. Error z value Pr(>|z|)    
## (Intercept)              0.7733     0.1101   7.024 2.15e-12 ***
## Invertebrate_predators   0.3146     0.2403   1.309    0.191    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(2.2127) family taken to be 1)
## 
##     Null deviance: 149.88  on 119  degrees of freedom
## Residual deviance: 148.05  on 118  degrees of freedom
## AIC: 488.88
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  2.213 
##           Std. Err.:  0.701 
## 
##  2 x log-likelihood:  -482.884

AICc(my_glm)

## [1] 489.0908

Abundance of fish predators

my_glm<- glm.nb(Consumed_individuals~ Predatory_fish, data=Squidpops2_stand)
summary(my_glm)

## 
## Call:
## glm.nb(formula = Consumed_individuals ~ Predatory_fish, data = Squidpops2_stand, 
##     init.theta = 2.124376703, link = log)
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      0.8129     0.1299   6.259 3.86e-10 ***
## Predatory_fish   0.1349     0.2519   0.535    0.592    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(2.1244) family taken to be 1)
## 
##     Null deviance: 147.65  on 119  degrees of freedom
## Residual deviance: 147.34  on 118  degrees of freedom
## AIC: 490.38
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  2.124 
##           Std. Err.:  0.655 
## 
##  2 x log-likelihood:  -484.377

AICc(my_glm)

## [1] 490.5841

Test interaction between temperature and dissolved oxygen

my_glmInt<- glm.nb(Consumed_individuals~ Temperature*D_oxygen_mgL, data=Squidpops2_stand)
summary(my_glmInt)

## 
## Call:
## glm.nb(formula = Consumed_individuals ~ Temperature * D_oxygen_mgL, 
##     data = Squidpops2_stand, init.theta = 3.311609876, link = log)
## 
## Coefficients:
##                          Estimate Std. Error z value Pr(>|z|)   
## (Intercept)                0.8111     0.2619   3.097  0.00196 **
## Temperature               -0.5330     0.3853  -1.383  0.16660   
## D_oxygen_mgL              -0.4962     0.7885  -0.629  0.52917   
## Temperature:D_oxygen_mgL   1.7523     0.8777   1.997  0.04588 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(3.3116) family taken to be 1)
## 
##     Null deviance: 171.19  on 119  degrees of freedom
## Residual deviance: 154.00  on 116  degrees of freedom
## AIC: 479.3
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  3.31 
##           Std. Err.:  1.36 
## 
##  2 x log-likelihood:  -469.298

AICc(my_glmInt)

## [1] 479.8243

5f) Model selection: Squidpops 24h

Temperature

my_glmTS24<- glm(Consumed_individuals~  Temperature, data=Squidpops24_stand, family = poisson)
check_overdispersion(my_glmTS24)

## # Overdispersion test
## 
##        dispersion ratio =  0.331
##   Pearson's Chi-Squared = 39.057
##                 p-value =      1

## No overdispersion detected.

summary(my_glmTS24)

## 
## Call:
## glm(formula = Consumed_individuals ~ Temperature, family = poisson, 
##     data = Squidpops24_stand)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   1.1227     0.1166   9.632  < 2e-16 ***
## Temperature   0.4879     0.1575   3.098  0.00195 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 61.719  on 119  degrees of freedom
## Residual deviance: 51.673  on 118  degrees of freedom
## AIC: 443.43
## 
## Number of Fisher Scoring iterations: 4

AICc(my_glmTS24)

## [1] 443.5292

Salinity

my_glmSS24<- glm(Consumed_individuals~  Salinity_ppt, data=Squidpops24_stand, family = poisson)
summary(my_glmSS24)

## 
## Call:
## glm(formula = Consumed_individuals ~ Salinity_ppt, family = poisson, 
##     data = Squidpops24_stand)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    1.1271     0.1434   7.857 3.93e-15 ***
## Salinity_ppt   0.4670     0.1956   2.387    0.017 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 61.719  on 119  degrees of freedom
## Residual deviance: 55.707  on 118  degrees of freedom
## AIC: 447.46
## 
## Number of Fisher Scoring iterations: 4

AICc(my_glmSS24)

## [1] 447.5637

Dissolved oxygen

my_glm<- glm(Consumed_individuals~  D_oxygen_mgL, data=Squidpops24_stand, family = poisson)
summary(my_glm)

## 
## Call:
## glm(formula = Consumed_individuals ~ D_oxygen_mgL, family = poisson, 
##     data = Squidpops24_stand)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   1.39170    0.08033  17.324   <2e-16 ***
## D_oxygen_mgL  0.12467    0.14949   0.834    0.404    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 61.719  on 119  degrees of freedom
## Residual deviance: 61.028  on 118  degrees of freedom
## AIC: 452.78
## 
## Number of Fisher Scoring iterations: 4

AICc(my_glm)

## [1] 452.8848

Visibility

my_glm<- glm(Consumed_individuals~ V_Visibility, data=Squidpops24_stand, family = poisson)
summary(my_glm)

## 
## Call:
## glm(formula = Consumed_individuals ~ V_Visibility, family = poisson, 
##     data = Squidpops24_stand)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   1.53164    0.08209  18.659   <2e-16 ***
## V_Visibility -0.20381    0.16916  -1.205    0.228    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 61.719  on 119  degrees of freedom
## Residual deviance: 60.248  on 118  degrees of freedom
## AIC: 452
## 
## Number of Fisher Scoring iterations: 4

AICc(my_glm)

## [1] 452.1049

Community composition

my_glm<- glm(Consumed_individuals~  PCo1, data=Squidpops24_stand, family = poisson)
summary(my_glm)

## 
## Call:
## glm(formula = Consumed_individuals ~ PCo1, family = poisson, 
##     data = Squidpops24_stand)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.50543    0.06171  24.394   <2e-16 ***
## PCo1        -0.16872    0.12776  -1.321    0.187    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 61.719  on 119  degrees of freedom
## Residual deviance: 59.944  on 118  degrees of freedom
## AIC: 451.7
## 
## Number of Fisher Scoring iterations: 4

AICc(my_glm)

## [1] 451.8008

Kelp cover

my_glmK<- glm(Consumed_individuals~ Kelp, data=Squidpops24_stand, family = poisson)
summary(my_glmK)

## 
## Call:
## glm(formula = Consumed_individuals ~ Kelp, family = poisson, 
##     data = Squidpops24_stand)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.37365    0.06285  21.856   <2e-16 ***
## Kelp         0.25829    0.15153   1.705   0.0883 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 61.719  on 119  degrees of freedom
## Residual deviance: 58.904  on 118  degrees of freedom
## AIC: 450.66
## 
## Number of Fisher Scoring iterations: 4

AICc(my_glmK)

## [1] 450.7603

Bare rock cover

my_glm<- glm(Consumed_individuals~ Bare_rock, data=Squidpops24_stand, family = poisson)
summary(my_glm)

## 
## Call:
## glm(formula = Consumed_individuals ~ Bare_rock, family = poisson, 
##     data = Squidpops24_stand)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.51985    0.06556  23.183   <2e-16 ***
## Bare_rock   -0.19684    0.13441  -1.464    0.143    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 61.719  on 119  degrees of freedom
## Residual deviance: 59.542  on 118  degrees of freedom
## AIC: 451.3
## 
## Number of Fisher Scoring iterations: 4

AICc(my_glm)

## [1] 451.3986

Species richness

my_glm<- glm(Consumed_individuals~ Mean_SR, data=Squidpops24_stand, family = poisson)
summary(my_glm)

## 
## Call:
## glm(formula = Consumed_individuals ~ Mean_SR, family = poisson, 
##     data = Squidpops24_stand)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.40627    0.08759  16.054   <2e-16 ***
## Mean_SR      0.08470    0.15635   0.542    0.588    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 61.719  on 119  degrees of freedom
## Residual deviance: 61.426  on 118  degrees of freedom
## AIC: 453.18
## 
## Number of Fisher Scoring iterations: 4

AICc(my_glm)

## [1] 453.2822

Abundance of invertebrate predators

my_glm<- glm(Consumed_individuals~ Invertebrate_predators, data=Squidpops24_stand, family = poisson)
summary(my_glm)

## 
## Call:
## glm(formula = Consumed_individuals ~ Invertebrate_predators, 
##     family = poisson, data = Squidpops24_stand)
## 
## Coefficients:
##                        Estimate Std. Error z value Pr(>|z|)    
## (Intercept)             1.41039    0.05698  24.754   <2e-16 ***
## Invertebrate_predators  0.13100    0.12495   1.048    0.294    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 61.719  on 119  degrees of freedom
## Residual deviance: 60.641  on 118  degrees of freedom
## AIC: 452.39
## 
## Number of Fisher Scoring iterations: 4

AICc(my_glm)

## [1] 452.4973

Abundance of predatory fish

my_glm<- glm(Consumed_individuals~ Predatory_fish, data=Squidpops24_stand, family = poisson)
summary(my_glm)

## 
## Call:
## glm(formula = Consumed_individuals ~ Predatory_fish, family = poisson, 
##     data = Squidpops24_stand)
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     1.37803    0.06749  20.418   <2e-16 ***
## Predatory_fish  0.17727    0.12762   1.389    0.165    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 61.719  on 119  degrees of freedom
## Residual deviance: 59.815  on 118  degrees of freedom
## AIC: 451.57
## 
## Number of Fisher Scoring iterations: 4

AICc(my_glm)

## [1] 451.6715

Fish length (cm)

my_glmL<- glm(Consumed_individuals~ Fish_size, data=Squidpops24_stand, family = poisson)
summary(my_glmL)

## 
## Call:
## glm(formula = Consumed_individuals ~ Fish_size, family = poisson, 
##     data = Squidpops24_stand)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   1.3041     0.1144  11.400   <2e-16 ***
## Fish_size     0.2283     0.1665   1.371     0.17    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 61.719  on 119  degrees of freedom
## Residual deviance: 59.809  on 118  degrees of freedom
## AIC: 451.56
## 
## Number of Fisher Scoring iterations: 4

AICc(my_glmL)

## [1] 451.6657

Test interaction between temperature and dissolved oxygen

my_glm<- glm(Consumed_individuals~ Temperature * D_oxygen_mgL, data=Squidpops24_stand, family = poisson)
summary(my_glm)

## 
## Call:
## glm(formula = Consumed_individuals ~ Temperature * D_oxygen_mgL, 
##     family = poisson, data = Squidpops24_stand)
## 
## Coefficients:
##                          Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                1.2170     0.1628   7.474 7.78e-14 ***
## Temperature                0.4416     0.2252   1.961   0.0499 *  
## D_oxygen_mgL              -0.4196     0.4782  -0.877   0.3803    
## Temperature:D_oxygen_mgL   0.3665     0.5390   0.680   0.4965    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 61.719  on 119  degrees of freedom
## Residual deviance: 50.734  on 116  degrees of freedom
## AIC: 446.49
## 
## Number of Fisher Scoring iterations: 4

AICc(my_glm)

## [1] 446.8355

Part 6: Herbivory intensity: Effect of site specific variables

First, create datasets to be used in the analysis

Herbivory_analysis <- Herbivory_6 %>%
  dplyr::select(-Rope, -Bait, -Camera, -Hours, -Initial_length_mm, -Final_Length_mm, -Scraped, -Survival_proportion)

1) Standardise

Herb_all_num <- Herbivory_analysis %>%
  dplyr::select(-Site, -Bait2,-Consumed_proportion, -Day_length) 
Herb_all_stand <- decostand(Herb_all_num, method="range")
#Test autocorrelation 
corrplot(cor(Herb_all_stand), method = "number")

Herb_all_stand$Site <- Herbivory_analysis$Site
Herb_all_stand$Consumption<- Herbivory_analysis$Consumption
Herb_all_stand$Consumed_proportion<- Herbivory_analysis$Consumed_proportion

6a) Model selection

Temperature

my_glmTL<- glm(Consumed_proportion ~  Temperature, data=Herb_all_stand, family = binomial)

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

summary(my_glmTL)

## 
## Call:
## glm(formula = Consumed_proportion ~ Temperature, family = binomial, 
##     data = Herb_all_stand)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.7420     0.2145  -3.459 0.000542 ***
## Temperature  -0.7899     0.3171  -2.491 0.012738 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 491.97  on 589  degrees of freedom
## Residual deviance: 485.87  on 588  degrees of freedom
## AIC: 641.16
## 
## Number of Fisher Scoring iterations: 4

AICc(my_glmTL)

## [1] 641.184

check_overdispersion(my_glmTL)

## # Overdispersion test
## 
##  dispersion ratio =   0.784
##           p-value = < 0.001

## Underdispersion detected.

Salinity

my_glmSL<- glm(Consumed_proportion ~  Salinity_ppt, data=Herb_all_stand, family = binomial)

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

summary(my_glmSL)

## 
## Call:
## glm(formula = Consumed_proportion ~ Salinity_ppt, family = binomial, 
##     data = Herb_all_stand)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    0.8452     0.2733   3.093  0.00198 ** 
## Salinity_ppt  -3.3000     0.4135  -7.980 1.46e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 491.97  on 589  degrees of freedom
## Residual deviance: 418.33  on 588  degrees of freedom
## AIC: 575.41
## 
## Number of Fisher Scoring iterations: 4

AICc(my_glmSL)

## [1] 575.4276

Dissolved oxygen

my_glmOxL<- glm(Consumed_proportion ~  D_oxygen_mgL, data=Herb_all_stand, family = binomial)

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

summary(my_glmOxL)

## 
## Call:
## glm(formula = Consumed_proportion ~ D_oxygen_mgL, family = binomial, 
##     data = Herb_all_stand)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -1.9163     0.1931  -9.922  < 2e-16 ***
## D_oxygen_mgL   1.4687     0.3319   4.425 9.63e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 491.97  on 589  degrees of freedom
## Residual deviance: 472.31  on 588  degrees of freedom
## AIC: 624.7
## 
## Number of Fisher Scoring iterations: 4

AICc(my_glmOxL)

## [1] 624.7162

Visibility

my_glm<- glm(Consumed_proportion ~  V_Visibility, data=Herb_all_stand, family = binomial)

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

summary(my_glm)

## 
## Call:
## glm(formula = Consumed_proportion ~ V_Visibility, family = binomial, 
##     data = Herb_all_stand)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.19110    0.18540  -6.424 1.32e-10 ***
## V_Visibility -0.09443    0.37397  -0.253    0.801    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 491.97  on 589  degrees of freedom
## Residual deviance: 491.91  on 588  degrees of freedom
## AIC: 648.94
## 
## Number of Fisher Scoring iterations: 4

AICc(my_glm)

## [1] 648.9651

Community composition

my_glm<- glm(Consumed_proportion ~  PCo1, data=Herb_all_stand, family = binomial)

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

summary(my_glm)

## 
## Call:
## glm(formula = Consumed_proportion ~ PCo1, family = binomial, 
##     data = Herb_all_stand)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.6529     0.1519 -10.885  < 2e-16 ***
## PCo1          1.0650     0.2652   4.016 5.91e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 491.97  on 589  degrees of freedom
## Residual deviance: 476.04  on 588  degrees of freedom
## AIC: 633.32
## 
## Number of Fisher Scoring iterations: 4

AICc(my_glm)

## [1] 633.3399

Kelp cover

my_glm<- glm(Consumed_proportion ~  Kelp, data=Herb_all_stand, family = binomial)

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

summary(my_glm)

## 
## Call:
## glm(formula = Consumed_proportion ~ Kelp, family = binomial, 
##     data = Herb_all_stand)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.9501     0.1366  -6.957 3.49e-12 ***
## Kelp         -1.1110     0.4085  -2.720  0.00653 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 491.97  on 589  degrees of freedom
## Residual deviance: 483.70  on 588  degrees of freedom
## AIC: 638.16
## 
## Number of Fisher Scoring iterations: 4

AICc(my_glm)

## [1] 638.1771

Bare rock cover

my_glmBRL<- glm(Consumed_proportion ~  Bare_rock, data=Herb_all_stand, family = binomial)

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

summary(my_glmBRL)

## 
## Call:
## glm(formula = Consumed_proportion ~ Bare_rock, family = binomial, 
##     data = Herb_all_stand)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -2.0942     0.1778 -11.780  < 2e-16 ***
## Bare_rock     1.9712     0.2940   6.704 2.03e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 491.97  on 589  degrees of freedom
## Residual deviance: 445.04  on 588  degrees of freedom
## AIC: 603.94
## 
## Number of Fisher Scoring iterations: 4

AICc(my_glmBRL)

## [1] 603.9645

Species richness

my_glm<- glm(Consumed_proportion ~  Mean_SR, data=Herb_all_stand, family = binomial)

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

summary(my_glm)

## 
## Call:
## glm(formula = Consumed_proportion ~ Mean_SR, family = binomial, 
##     data = Herb_all_stand)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.3117     0.1946  -6.740 1.58e-11 ***
## Mean_SR       0.1688     0.3484   0.485    0.628    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 491.97  on 589  degrees of freedom
## Residual deviance: 491.74  on 588  degrees of freedom
## AIC: 647.32
## 
## Number of Fisher Scoring iterations: 4

AICc(my_glm)

## [1] 647.339

Abundance of herbivores

my_glm<- glm(Consumed_proportion ~  Herbivores, data=Herb_all_stand, family = binomial)

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

summary(my_glm)

## 
## Call:
## glm(formula = Consumed_proportion ~ Herbivores, family = binomial, 
##     data = Herb_all_stand)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.2632     0.1225 -10.308   <2e-16 ***
## Herbivores    0.1492     0.3326   0.449    0.654    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 491.97  on 589  degrees of freedom
## Residual deviance: 491.77  on 588  degrees of freedom
## AIC: 647.94
## 
## Number of Fisher Scoring iterations: 4

AICc(my_glm)

## [1] 647.9606

Phlorotannins

my_glmPH<- glm(Consumed_proportion ~  phloro_mean, data=Herb_all_stand, family = binomial)

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

summary(my_glmPH)

## 
## Call:
## glm(formula = Consumed_proportion ~ phloro_mean, family = binomial, 
##     data = Herb_all_stand)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.2487     0.1392  -1.787   0.0739 .  
## phloro_mean  -4.5287     0.6353  -7.128 1.02e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 491.97  on 589  degrees of freedom
## Residual deviance: 398.55  on 588  degrees of freedom
## AIC: 544.61
## 
## Number of Fisher Scoring iterations: 6

AICc(my_glmPH)

## [1] 544.6339

Part 7: Visualisation of best fitted predictors

7a): Crab consumption

Make predictions according to individual models with lowest AICc

Crabs 2 h

newdfrC2hT = expand.grid(Temperature = rep(seq(from = min(Crabs2_stand$Temperature), to = max(Crabs2_stand$Temperature), length.out = 100), 5))
 
#add preditions
dat1 <- predict(my_glmT, newdata = newdfrC2hT, type = "link", se.fit = TRUE)
critval <- 1.96 ## approx 95% CI
upr <- dat1$fit + (critval * dat1$se.fit)
lwr <- dat1$fit - (critval * dat1$se.fit)
fit <- dat1$fit
fit2 <- my_glmT$family$linkinv(fit)
upr2 <- my_glmT$family$linkinv(upr)
lwr2 <- my_glmT$family$linkinv(lwr)

newdfrC2hT$lwr <- lwr2 
newdfrC2hT$upr <- upr2 
newdfrC2hT$fit <- fit2
#For Salinity

newdfrC2S = expand.grid(D_oxygen_mgL = rep(seq(from = min(Crabs2_stand$D_oxygen_mgL), to = max(Crabs2_stand$D_oxygen_mgL), length.out = 100), 5))
 
dat2 <- predict(my_glmOx, newdata = newdfrC2S, type = "link", se.fit = TRUE)
upr1 <- dat2$fit + (critval * dat2$se.fit)
lwr1 <- dat2$fit - (critval * dat2$se.fit)
fit1 <- dat2$fit
fit3 <- my_glmOx$family$linkinv(fit1)
upr3 <- my_glmOx$family$linkinv(upr1)
lwr3 <- my_glmOx$family$linkinv(lwr1)

newdfrC2S$lwr <- lwr3 
newdfrC2S$upr <- upr3 
newdfrC2S$fit <- fit3


#For Bare rock cover 

newdfrC2BR = expand.grid(Bare_rock = rep(seq(from = min(Crabs2_stand$Bare_rock), to = max(Crabs2_stand$Bare_rock), length.out = 100), 5))
 
dat2 <- predict(my_glmBR, newdata = newdfrC2BR , type = "link", se.fit = TRUE)
upr1 <- dat2$fit + (critval * dat2$se.fit)
lwr1 <- dat2$fit - (critval * dat2$se.fit)
fit1 <- dat2$fit
fit3 <- my_glmBR$family$linkinv(fit1)
upr3 <- my_glmBR$family$linkinv(upr1)
lwr3 <- my_glmBR$family$linkinv(lwr1)

newdfrC2BR$lwr <- lwr3 
newdfrC2BR$upr <- upr3 
newdfrC2BR$fit <- fit3

#For Salinity

newdfrC2Sal = expand.grid(Salinity_ppt = rep(seq(from = min(Crabs2_stand$Salinity_ppt), to = max(Crabs2_stand$Salinity_ppt), length.out = 100), 5))
 
dat2 <- predict(my_glmSC2, newdata = newdfrC2Sal , type = "link", se.fit = TRUE)
upr1 <- dat2$fit + (critval * dat2$se.fit)
lwr1 <- dat2$fit - (critval * dat2$se.fit)
fit1 <- dat2$fit
fit3 <- my_glmSC2$family$linkinv(fit1)
upr3 <- my_glmSC2$family$linkinv(upr1)
lwr3 <- my_glmSC2$family$linkinv(lwr1)

newdfrC2Sal$lwr <- lwr3 
newdfrC2Sal$upr <- upr3 
newdfrC2Sal$fit <- fit3

Crabs 24 h

newdfrC24T = expand.grid(Temperature = rep(seq(from = min(Crabs24_stand$Temperature), to = max(Crabs24_stand$Temperature), length.out = 100), 5))
 
#add preditions
dat1 <- predict(my_glmT24, newdata = newdfrC24T, type = "link", se.fit = TRUE)
critval <- 1.96 ## approx 95% CI
upr <- dat1$fit + (critval * dat1$se.fit)
lwr <- dat1$fit - (critval * dat1$se.fit)
fit <- dat1$fit
fit2 <- my_glmT24$family$linkinv(fit)
upr2 <- my_glmT24$family$linkinv(upr)
lwr2 <- my_glmT24$family$linkinv(lwr)

newdfrC24T$lwr <- lwr2 
newdfrC24T$upr <- upr2 
newdfrC24T$fit <- fit2
#For Salinity

newdfrC24S = expand.grid(Salinity_ppt = rep(seq(from = min(Crabs24_stand$Salinity_ppt), to = max(Crabs24_stand$Salinity_ppt), length.out = 100), 5))
 
dat2 <- predict(my_glmS24, newdata = newdfrC24S, type = "link", se.fit = TRUE)
upr1 <- dat2$fit + (critval * dat2$se.fit)
lwr1 <- dat2$fit - (critval * dat2$se.fit)
fit1 <- dat2$fit
fit3 <- my_glmS24$family$linkinv(fit1)
upr3 <- my_glmS24$family$linkinv(upr1)
lwr3 <- my_glmS24$family$linkinv(lwr1)

newdfrC24S$lwr <- lwr3 
newdfrC24S$upr <- upr3 
newdfrC24S$fit <- fit3


#For Bare rock cover 

newdfrC24BR = expand.grid(Bare_rock = rep(seq(from = min(Crabs24_stand$Bare_rock), to = max(Crabs24_stand$Bare_rock), length.out = 100), 5))
 
dat2 <- predict(my_glmBR24, newdata = newdfrC24BR , type = "link", se.fit = TRUE)
upr1 <- dat2$fit + (critval * dat2$se.fit)
lwr1 <- dat2$fit - (critval * dat2$se.fit)
fit1 <- dat2$fit
fit3 <- my_glmBR24$family$linkinv(fit1)
upr3 <- my_glmBR24$family$linkinv(upr1)
lwr3 <- my_glmBR24$family$linkinv(lwr1)

newdfrC24BR$lwr <- lwr3 
newdfrC24BR$upr <- upr3 
newdfrC24BR$fit <- fit3

#For Fish length 

newdfrC24FL = expand.grid(Fish_size = rep(seq(from = min(Crabs24_stand$Fish_size), to = max(Crabs24_stand$Fish_size), length.out = 100), 5))
 
dat2 <- predict(my_glm24_FL, newdata = newdfrC24FL , type = "link", se.fit = TRUE)
upr1 <- dat2$fit + (critval * dat2$se.fit)
lwr1 <- dat2$fit - (critval * dat2$se.fit)
fit1 <- dat2$fit
fit3 <- my_glm24_FL$family$linkinv(fit1)
upr3 <- my_glm24_FL$family$linkinv(upr1)
lwr3 <- my_glm24_FL$family$linkinv(lwr1)

newdfrC24FL$lwr <- lwr3 
newdfrC24FL$upr <- upr3 
newdfrC24FL$fit <- fit3

1) Temperature vs Crabs

ggplot() +
  geom_ribbon(data = newdfrC24T, aes(x = Temperature, y = fit, ymin = lwr, ymax = upr), fill = "orange1", alpha = .25) +
  geom_line(data = newdfrC24T, aes(x = Temperature, y =fit), color =  "orange1", linewidth = 0.6) +
  geom_jitter(data = Crabs24_stand, aes(x= Temperature, y=Consumed_individuals), color = "orange1", width=0.01, height = 0.1, size = 1, alpha = 0.5)+
  geom_ribbon(data = newdfrC2hT, aes(x = Temperature, y = fit, ymin = lwr, ymax = upr), fill = "darkslateblue", alpha = .25) +
  geom_line(data = newdfrC2hT, aes(x = Temperature, y =fit), color =  "darkslateblue", linewidth = 0.6) +
  geom_jitter(data = Crabs2_stand, aes(x= Temperature, y=Consumed_individuals), color = "darkslateblue", width=0.01, height = 0.1, size = 1, alpha = 0.5)+
  labs(x = "Temperature (°C)", y = "Crabs consumed")  + Theme1 +
  scale_x_continuous(breaks = c(-0.04868, 0.21096, 0.4706, 0.73023, 0.98987), labels = c("10", "12", "14", "16", "18"), limits = c(0,1)) +
  scale_y_continuous(limits = c(0,6.3), breaks = seq(0, 5, 1))

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

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

ggsave("TempC.jpeg", units="mm", width=110, height=80, dpi=300)

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

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

2) Salinity vs Crabs

ggplot() +
  geom_ribbon(data = newdfrC24S, aes(x = Salinity_ppt, y = fit, ymin = lwr, ymax = upr), fill = "orange1", alpha = .25) +
  geom_line(data = newdfrC24S, aes(x = Salinity_ppt, y =fit), color =  "orange1", size = 1) +
  geom_jitter(data = Crabs24_stand, aes(x= Salinity_ppt, y=Consumed_individuals), color = "orange1", width=0.01, height = 0.1, size = 1.5, alpha = 0.5)+
  geom_ribbon(data = newdfrC2Sal, aes(x = Salinity_ppt, y = fit, ymin = lwr, ymax = upr), fill = "darkslateblue", alpha = .25) +
  geom_line(data = newdfrC2Sal, aes(x = Salinity_ppt, y =fit), color =  "darkslateblue", size = 1) +
  geom_jitter(data = Crabs2_stand, aes(x= Salinity_ppt, y=Consumed_individuals), color = "darkslateblue", width=0.01, height = 0.1, size = 1.5, alpha = 0.5)+
  labs(x = "Salinity (ppt)", y = "Crabs consumed")  + Theme1 +
  scale_x_continuous(breaks = c(0.09803922, 0.49019608, 0.88235299, 1.2745098), labels = c("33", "33.5", "34", "33.5"), limits = c(0,1)) +
  scale_y_continuous(limits = c(0,5.8), breaks = seq(0, 5, 1)) + 
  theme(legend.position = c(0.1, 0.85))

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

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

ggsave("SalC.jpeg", units="mm", width=110, height=90, dpi=300)

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

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

3) Oxygen vs Crabs

ggplot() +
  geom_ribbon(data = newdfrC2S, aes(x = D_oxygen_mgL, y = fit, ymin = lwr, ymax = upr), fill = "darkslateblue", alpha = .25) +
  geom_line(data = newdfrC2S, aes(x = D_oxygen_mgL, y =fit), color =  "darkslateblue", size = 1) +
  geom_jitter(data = Crabs2_stand, aes(x= D_oxygen_mgL, y=Consumed_individuals), color = "darkslateblue", width=0.01, height = 0.1, size = 1.5, alpha = 0.5)+
  geom_jitter(data = Crabs24_stand, aes(x= D_oxygen_mgL, y=Consumed_individuals), color = "orange1", width=0.01, height = 0.1, size = 1.5, alpha = 0.5)+
  labs(x = "Dissolved oxygen (mgL)", y = "Crabs consumed")  + Theme1 +
  scale_x_continuous(breaks = c(0.1082296, 0.384854, 0.66147994, 0.935517247), labels = c("4", "6", "8", "10"), limits = c(0,1)) +
  scale_y_continuous(limits = c(0,5.8), breaks = seq(0, 5, 1)) + 
  theme(legend.position = c(0.1, 0.85))

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

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

ggsave("OxC.jpeg", units="mm", width=110, height=80, dpi=300)

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

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

4) Bare rock Crabs 24 h

ggplot() +
  geom_ribbon(data = newdfrC24BR, aes(x = Bare_rock, y = fit, ymin = lwr, ymax = upr), fill = "orange1", alpha = .25) +
  geom_line(data = newdfrC24BR, aes(x = Bare_rock, y =fit), color =  "orange1", size = 1) +
  geom_jitter(data = Crabs24_stand, aes(x= Bare_rock, y=Consumed_individuals), color = "orange1", width=0.02, height = 0.1, size = 1.5, alpha = 0.5)+
  geom_ribbon(data = newdfrC2BR, aes(x = Bare_rock, y = fit, ymin = lwr, ymax = upr), fill = "darkslateblue", alpha = .25) +
  geom_line(data = newdfrC2BR, aes(x = Bare_rock, y =fit), color =  "darkslateblue", size = 1) +
  geom_jitter(data = Crabs2_stand, aes(x= Bare_rock, y=Consumed_individuals), color = "darkslateblue", width=0.02, height = 0.1, size = 1.5, alpha = 0.5)+
  labs(x = "Bare rock cover (%)", y = "Consumed crabs")  + Theme1 +
  scale_x_continuous(breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1), labels = c("0", "20", "40", "60", "80", "100"), limits = c(0,1)) +
  scale_y_continuous(limits = c(0,5.8), breaks = seq(0, 5, 1)) + 
  theme(legend.position = c(0.1, 0.85))

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

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

ggsave("BareC.jpeg", units="mm", width=110, height=80, dpi=300)

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

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

4) Fish length - Crabs 24 h

ggplot() +
  geom_ribbon(data = newdfrC24FL, aes(x = Fish_size, y = fit, ymin = lwr, ymax = upr), fill = "orange1", alpha = .25) +
  geom_line(data = newdfrC24FL, aes(x = Fish_size, y =fit), color =  "orange1", size = 1) +
  geom_jitter(data = Crabs24_stand, aes(x= Fish_size, y=Consumed_individuals), color = "orange1", width=0.02, height = 0.1, size = 1.5, alpha = 0.5)+
  geom_jitter(data = Crabs2_stand, aes(x= Bare_rock, y=Consumed_individuals), color = "darkslateblue", width=0.02, height = 0.1, size = 1.5, alpha = 0.5)+
  labs(x = "Fish length (cm)", y = "Consumed crabs")  + Theme1 +
  scale_x_continuous(breaks = c(0.011979959
,0.165667696
, 0.319355434
, 0.473043171
, 0.626730908
, 0.780418645
,0.934106383), labels = c("14","16", "18", "20", "22", "24", "26"), limits = c(0,1)) +
  scale_y_continuous(limits = c(0,5.8), breaks = seq(0, 5, 1)) + 
  theme(legend.position = c(0.1, 0.85))

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

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

ggsave("Fish_length.jpeg", units="mm", width=110, height=80, dpi=300)

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

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

7b): Squidpop consumption

Squidpops 2 h

S2Ox = expand.grid(D_oxygen_mgL = rep(seq(from = min(Squidpops2_stand$D_oxygen_mgL), to = max(Squidpops2_stand$D_oxygen_mgL), length.out = 100), 5))
 
#add preditions
dat1 <- predict(my_glmOxS, newdata = S2Ox, type = "link", se.fit = TRUE)
critval <- 1.96 ## approx 95% CI
upr <- dat1$fit + (critval * dat1$se.fit)
lwr <- dat1$fit - (critval * dat1$se.fit)
fit <- dat1$fit
fit2 <- my_glmOxS$family$linkinv(fit)
upr2 <- my_glmOxS$family$linkinv(upr)
lwr2 <- my_glmOxS$family$linkinv(lwr)

S2Ox$lwr <- lwr2 
S2Ox$upr <- upr2 
S2Ox$fit <- fit2
#For Pco

S2PCo1 = expand.grid(PCo1 = rep(seq(from = min(Squidpops2_stand$PCo1), to = max(Squidpops2_stand$PCo1), length.out = 100), 5))
 
#add preditions
dat1 <- predict(my_glmPC, newdata = S2PCo1, type = "link", se.fit = TRUE)
critval <- 1.96 ## approx 95% CI
upr <- dat1$fit + (critval * dat1$se.fit)
lwr <- dat1$fit - (critval * dat1$se.fit)
fit <- dat1$fit
fit2 <- my_glmPC$family$linkinv(fit)
upr2 <- my_glmPC$family$linkinv(upr)
lwr2 <- my_glmPC$family$linkinv(lwr)

S2PCo1$lwr <- lwr2 
S2PCo1$upr <- upr2 
S2PCo1$fit <- fit2


#For Fish length


S2FL = expand.grid(Fish_size = rep(seq(from = min(Squidpops2_stand$Fish_size), to = max(Squidpops2_stand$Fish_size), length.out = 100), 5))
 
#add preditions
dat1 <- predict(my_glmL, newdata = S2FL, type = "link", se.fit = TRUE)
critval <- 1.96 ## approx 95% CI
upr <- dat1$fit + (critval * dat1$se.fit)
lwr <- dat1$fit - (critval * dat1$se.fit)
fit <- dat1$fit
fit2 <- my_glmL$family$linkinv(fit)
upr2 <- my_glmL$family$linkinv(upr)
lwr2 <- my_glmL$family$linkinv(lwr)

S2FL$lwr <- lwr2 
S2FL$upr <- upr2 
S2FL$fit <- fit2

Squidpops 24 h

S24T = expand.grid(Temperature = rep(seq(from = min(Squidpops24_stand$Temperature), to = max(Squidpops24_stand$Temperature), length.out = 100), 5))
 
#add preditions
dat1 <- predict(my_glmTS24, newdata = S24T, type = "link", se.fit = TRUE)
critval <- 1.96 ## approx 95% CI
upr <- dat1$fit + (critval * dat1$se.fit)
lwr <- dat1$fit - (critval * dat1$se.fit)
fit <- dat1$fit
fit2 <- my_glmTS24$family$linkinv(fit)
upr2 <- my_glmTS24$family$linkinv(upr)
lwr2 <- my_glmTS24$family$linkinv(lwr)

S24T$lwr <- lwr2 
S24T$upr <- upr2 
S24T$fit <- fit2
#For Salinity

S24Sal = expand.grid(Salinity_ppt = rep(seq(from = min(Squidpops24_stand$Salinity_ppt), to = max(Squidpops24_stand$Salinity_ppt), length.out = 100), 5))
 
#add preditions
dat1 <- predict(my_glmSS24, newdata = S24Sal, type = "link", se.fit = TRUE)
critval <- 1.96 ## approx 95% CI
upr <- dat1$fit + (critval * dat1$se.fit)
lwr <- dat1$fit - (critval * dat1$se.fit)
fit <- dat1$fit
fit2 <- my_glmSS24$family$linkinv(fit)
upr2 <- my_glmSS24$family$linkinv(upr)
lwr2 <- my_glmSS24$family$linkinv(lwr)

S24Sal$lwr <- lwr2 
S24Sal$upr <- upr2 
S24Sal$fit <- fit2

1) Temperature vs squidpops

ggplot() +
  geom_ribbon(data = S24T, aes(x = Temperature, y = fit, ymin = lwr, ymax = upr), fill = "orange1", alpha = .25) +
  geom_line(data = S24T, aes(x = Temperature, y =fit), color =  "orange1", size = 1) +
  geom_jitter(data = Squidpops24_stand, aes(x= Temperature, y=Consumed_individuals), color = "orange1", width=0.02, height = 0.1, size = 1.5, alpha = 0.5)+
   geom_jitter(data = Squidpops2_stand, aes(x= Temperature, y=Consumed_individuals), color = "darkslateblue", width=0.02, height = 0.1, size = 1.5, alpha = 0.5)+
  labs(x = "Temperature (°C)", y = "Consumed squidpops")  + Theme1 +
  scale_x_continuous(breaks = c(-0.04868, 0.21096, 0.4706, 0.73023, 0.98987), labels = c("10", "12", "14", "16", "18"), limits = c(0,1)) +
  scale_y_continuous(limits = c(0,5.8), breaks = seq(0, 5, 1)) + 
  theme(legend.position = c(0.1, 0.85))

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

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

ggsave("TemperatureS.jpeg", units="mm", width=110, height=80, dpi=300)

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

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

2) Salinity vs squidpops

ggplot() +
  geom_ribbon(data = S24Sal, aes(x = Salinity_ppt, y = fit, ymin = lwr, ymax = upr), fill = "orange1", alpha = .25) +
  geom_line(data = S24Sal, aes(x = Salinity_ppt, y =fit), color =  "orange1", size = 1) +
  geom_jitter(data = Squidpops24_stand, aes(x= Salinity_ppt, y=Consumed_individuals), color = "orange1", width=0.01, height = 0.1, size = 1.5, alpha = 0.5)+
   geom_jitter(data = Squidpops2_stand, aes(x= Salinity_ppt, y=Consumed_individuals), color = "darkslateblue", width=0.01, height = 0.1, size = 1.5, alpha = 0.5)+
  labs(x = "Salinity (ppt)", y = "Consumed squidpops")  + Theme1 +
  scale_x_continuous(breaks = c(0.09803922, 0.49019608, 0.88235299, 1.2745098), labels = c("33", "33.5", "34", "33.5"), limits = c(0,1)) +
  scale_y_continuous(limits = c(0,5.8), breaks = seq(0, 5, 1)) + 
  theme(legend.position = c(0.1, 0.85))

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

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

ggsave("SalS.jpeg", units="mm", width=110, height=80, dpi=300)

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

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

3) Oxygen vs squidpops

ggplot() +
  geom_ribbon(data = S2Ox, aes(x = D_oxygen_mgL, y = fit, ymin = lwr, ymax = upr), fill = "darkslateblue", alpha = .25) +
  geom_line(data = S2Ox, aes(x = D_oxygen_mgL, y =fit), color =  "darkslateblue", size = 1) +
  geom_jitter(data = Squidpops2_stand, aes(x= D_oxygen_mgL, y=Consumed_individuals), color = "darkslateblue", width=0.01, height = 0.1, size = 1.5, alpha = 0.5)+
  geom_jitter(data = Squidpops24_stand, aes(x= D_oxygen_mgL, y=Consumed_individuals), color = "orange1", width=0.01, height = 0.1, size = 1.5, alpha = 0.5)+
  labs(x = "Dissolved oxygen (mgL)", y = "Consumed squidpops")  + Theme1 +
  scale_x_continuous(breaks = c(0.1082296, 0.384854, 0.66147994, 0.935517247), labels = c("4", "6", "8", "10"), limits = c(0,1)) +
  scale_y_continuous(limits = c(0,5.8), breaks = seq(0, 5, 1)) + 
  theme(legend.position = c(0.1, 0.85))

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

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

ggsave("OxS.jpeg", units="mm", width=110, height=80, dpi=300)

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

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

Bare rock vs squidpops

ggplot() +
  geom_jitter(data = Squidpops2_stand, aes(x= Fish_size, y=Consumed_individuals), color = "darkslateblue", width=0.015, height = 0.1, size = 1.5, alpha = 0.5)+
    geom_jitter(data = Squidpops24_stand, aes(x= Fish_size, y=Consumed_individuals), color = "orange1", width=0.015, height = 0.1, size = 1.5, alpha = 0.5)+
  labs(x = "Bare rock cover (%)", y = "Consumed squidpops")  + Theme1 +
  scale_x_continuous(breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1), labels = c("0", "20", "40", "60", "80", "100"), limits = c(0,1)) +
  scale_y_continuous(limits = c(0,5.8), breaks = seq(0, 5, 1)) + 
  theme(legend.position = c(0.1, 0.85))

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

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

ggsave("BRS.jpeg", units="mm", width=110, height=80, dpi=300)

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

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

Fish length vs squidpops

ggplot() +
  geom_ribbon(data = S2FL, aes(x = Fish_size, y = fit, ymin = lwr, ymax = upr), fill = "darkslateblue", alpha = .25) +
  geom_line(data = S2FL, aes(x = Fish_size, y =fit), color =  "darkslateblue", size = 1) +
  geom_jitter(data = Squidpops2_stand, aes(x= Fish_size, y=Consumed_individuals), color = "darkslateblue", width=0.015, height = 0.1, size = 1.5, alpha = 0.5)+
    geom_jitter(data = Squidpops24_stand, aes(x= Fish_size, y=Consumed_individuals), color = "orange1", width=0.015, height = 0.1, size = 1.5, alpha = 0.5)+
  labs(x = "Fish Length (cm)", y = "Consumed squidpops")  + Theme1 +
  scale_x_continuous(breaks = c(0.011979959
,0.165667696
, 0.319355434
, 0.473043171
, 0.626730908
, 0.780418645
,0.934106383), labels = c("14","16", "18", "20", "22", "24", "26"), limits = c(0,1)) +
  scale_y_continuous(limits = c(0,5.8), breaks = seq(0, 5, 1)) + 
  theme(legend.position = c(0.1, 0.85))

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

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

ggsave("FLS2.jpeg", units="mm", width=110, height=80, dpi=300)

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

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

7c): Seaweed consumption

LPhlo = expand.grid(phloro_mean = rep(seq(from = min(Herb_all_stand$phloro_mean), to = max(Herb_all_stand$phloro_mean), length.out = 100), 5))
 
#add preditions
dat1 <- predict(my_glmPH, newdata = LPhlo, type = "link", se.fit = TRUE)
critval <- 1.96 ## approx 95% CI
upr <- dat1$fit + (critval * dat1$se.fit)
lwr <- dat1$fit - (critval * dat1$se.fit)
fit <- dat1$fit
fit2 <- my_glmPH$family$linkinv(fit)
upr2 <- my_glmPH$family$linkinv(upr)
lwr2 <- my_glmPH$family$linkinv(lwr)

LPhlo$lwr <- lwr2 
LPhlo$upr <- upr2 
LPhlo$fit <- fit2

#For Salinity
LSal = expand.grid(Salinity_ppt = rep(seq(from = min(Herb_all_stand$Salinity_ppt), to = max(Herb_all_stand$Salinity_ppt), length.out = 100), 5))
 
#add preditions
dat1 <- predict(my_glmSL, newdata = LSal, type = "link", se.fit = TRUE)
critval <- 1.96 ## approx 95% CI
upr <- dat1$fit + (critval * dat1$se.fit)
lwr <- dat1$fit - (critval * dat1$se.fit)
fit <- dat1$fit
fit2 <- my_glmSL$family$linkinv(fit)
upr2 <- my_glmSL$family$linkinv(upr)
lwr2 <- my_glmSL$family$linkinv(lwr)

LSal$lwr <- lwr2 
LSal$upr <- upr2 
LSal$fit <- fit2

#For Bare rock cover 

LBR = expand.grid(Bare_rock = rep(seq(from = min(Herb_all_stand$Bare_rock), to = max(Herb_all_stand$Bare_rock), length.out = 100), 5))
 
#add preditions
dat1 <- predict(my_glmBRL, newdata = LBR, type = "link", se.fit = TRUE)
critval <- 1.96 ## approx 95% CI
upr <- dat1$fit + (critval * dat1$se.fit)
lwr <- dat1$fit - (critval * dat1$se.fit)
fit <- dat1$fit
fit2 <- my_glmBRL$family$linkinv(fit)
upr2 <- my_glmBRL$family$linkinv(upr)
lwr2 <- my_glmBRL$family$linkinv(lwr)

LBR$lwr <- lwr2 
LBR$upr <- upr2 
LBR$fit <- fit2

#For Temperature

LTEMP = expand.grid(Temperature = rep(seq(from = min(Herb_all_stand$Temperature), to = max(Herb_all_stand$Temperature), length.out = 100), 5))
 
#add preditions
dat1 <- predict(my_glmTL, newdata = LTEMP, type = "link", se.fit = TRUE)
critval <- 1.96 ## approx 95% CI
upr <- dat1$fit + (critval * dat1$se.fit)
lwr <- dat1$fit - (critval * dat1$se.fit)
fit <- dat1$fit
fit2 <- my_glmTL$family$linkinv(fit)
upr2 <- my_glmTL$family$linkinv(upr)
lwr2 <- my_glmTL$family$linkinv(lwr)

LTEMP$lwr <- lwr2 
LTEMP$upr <- upr2 
LTEMP$fit <- fit2

#For oxygen

LOx = expand.grid(D_oxygen_mgL = rep(seq(from = min(Herb_all_stand$D_oxygen_mgL), to = max(Herb_all_stand$D_oxygen_mgL), length.out = 100), 5))
 
#add preditions
dat1 <- predict(my_glmOxL, newdata = LOx, type = "link", se.fit = TRUE)
critval <- 1.96 ## approx 95% CI
upr <- dat1$fit + (critval * dat1$se.fit)
lwr <- dat1$fit - (critval * dat1$se.fit)
fit <- dat1$fit
fit2 <- my_glmOxL$family$linkinv(fit)
upr2 <- my_glmOxL$family$linkinv(upr)
lwr2 <- my_glmOxL$family$linkinv(lwr)

LOx$lwr <- lwr2 
LOx$upr <- upr2 
LOx$fit <- fit2

Phlorotannins vs herbivory

ggplot() +
  geom_ribbon(data = LPhlo, aes(x = phloro_mean, y = fit, ymin = lwr, ymax = upr), fill = "orange1", alpha = .25) +
  geom_line(data = LPhlo, aes(x = phloro_mean, y =fit), color =  "orange1", size = 1) +
  geom_jitter(data = Herb_all_stand, aes(x= phloro_mean, y=Consumed_proportion), color = "orange1", width=0.03, height = 0.03, size = 1.5, alpha = 0.5)+
  labs(x = "Phlorotannin content (mg"~g^-1*" DW)", y = "Blade consumed proportion")  + Theme1 +
  scale_x_continuous(breaks = c(0.04788192, 0.2840443, 0.5202066, 0.756369, 0.9925314), labels = c("8", "16", "24", "32", "40"), limits = c(0,1)) +
  scale_y_continuous(limits = c(0,1.1), breaks = seq(0, 1, 0.2)) + 
  theme(legend.position = c(0.1, 0.85))

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

ggsave("PhloroL.jpeg", units="mm", width=110, height=80, dpi=300)

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

Salinity vs herbivory

ggplot() +
  geom_ribbon(data = LSal, aes(x = Salinity_ppt, y = fit, ymin = lwr, ymax = upr), fill = "orange1", alpha = .25) +
  geom_line(data = LSal, aes(x = Salinity_ppt, y =fit), color =  "orange1", size = 1) +
  geom_jitter(data = Herb_all_stand, aes(x= Salinity_ppt, y=Consumed_proportion), color = "orange1", width=0.03, height = 0.03, size = 1.5, alpha = 0.5)+
  labs(x = "Salinity (ppt)", y = "Blade consumed proportion")  + Theme1 +
  scale_x_continuous(breaks = c(0.09803922, 0.49019608, 0.88235299, 1.2745098), labels = c("33", "33.5", "34", "33.5"), limits = c(0,1)) +
  scale_y_continuous(limits = c(0,1.1), breaks = seq(0, 1, 0.2)) + 
  theme(legend.position = c(0.1, 0.85))

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

ggsave("SalL.jpeg", units="mm", width=110, height=80, dpi=300)

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

Bare rock vs herbivory

ggplot() +
  geom_ribbon(data = LBR, aes(x = Bare_rock, y = fit, ymin = lwr, ymax = upr), fill = "orange1", alpha = .25) +
  geom_line(data = LBR, aes(x = Bare_rock, y =fit), color =  "orange1", size = 1) +
  geom_jitter(data = Herb_all_stand, aes(x= Bare_rock, y=Consumed_proportion), color = "orange1", width=0.03, height = 0.03, size = 1.5, alpha = 0.5)+
  labs(x = "Bare rock cover (%)", y = "Predicted consumed proportion ")  + Theme1 +
  scale_x_continuous(breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1), labels = c("0", "20", "40", "60", "80", "100"), limits = c(0,1)) +
  scale_y_continuous(limits = c(0,1.1), breaks = seq(0, 1, 0.2)) + 
  theme(legend.position = c(0.1, 0.85))

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

ggsave("BareL.jpeg", units="mm", width=110, height=80, dpi=300)

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

Dissolved oxygen vs herbivory

ggplot() +
  geom_ribbon(data = LOx, aes(x = D_oxygen_mgL, y = fit, ymin = lwr, ymax = upr), fill = "orange1", alpha = .25) +
  geom_line(data = LOx, aes(x = D_oxygen_mgL, y =fit), color =  "orange1", linewidth = 0.6) +
  geom_jitter(data = Herb_all_stand, aes(x= D_oxygen_mgL, y=Consumed_proportion), color = "orange1", width=0.01, height = 0.03, size = 1.3, alpha = 0.5)+
  labs(x = "Dissolved oxygen (mgL)", y = "Predicted n° of consumed squidpops")  + Theme1 +
  scale_x_continuous(breaks = c(0.1082296, 0.384854, 0.66147994, 0.935517247), labels = c("4", "6", "8", "10"), limits = c(0,1))  +
  scale_y_continuous(limits = c(0,1.1), breaks = seq(0, 1, 0.2)) + 
  theme(legend.position = c(0.1, 0.85))

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

ggsave("DoxL.jpeg", units="mm", width=110, height=80, dpi=300)

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

Part 9: Supplementary tables

9a) Size estimations summary table

Size_summary <- Size_metadata %>%
  group_by(Species) %>%
  summarise (N_individuals = n(), mean_Size_sp = mean(Size_cm), sd = sd(Size_cm), Max_size = max(Size_cm), Min_size = min(Size_cm))

kable(Size_summary)

Species	N_individuals	mean_Size_sp	sd	Max_size	Min_size
Acanthistius_pictus	3	16.000000	1.0000000	17	15
Anisotremus_scapularis	4	21.250000	2.9860788	25	18
Aplodactylus_punctatus	54	27.944444	6.5196848	40	14
Auchenionchus_critinus	1	12.000000	NA	12	12
Auchenionchus_sp.	14	19.071429	6.0949448	32	11
Auchenionchus_variolosus	3	28.666667	5.6862407	35	24
Bovichtus_chilensis	3	17.666667	9.0184995	27	9
Calliclinus_geniguttatus	19	10.157895	2.3632172	14	4
Cheilodactylus_variegatus	125	22.560000	7.0715240	38	6
Cheilotrema_fasciatum	1	28.000000	NA	28	28
Chromis_crusma	85	13.611765	4.2457096	26	5
Congiopodus_peruvianus	2	13.500000	0.7071068	14	13
Girella_laevifrons	11	26.272727	7.3087743	38	15
Graus_nigra	15	25.066667	4.9923751	37	19
Hemilutjanus_macrophthalmos	2	20.000000	2.8284271	22	18
Hepatus_chilensis	1	6.000000	NA	6	6
Homalaspis_plana	20	9.800000	1.9084301	12	4
Hypsoblennius_sordidus	5	7.800000	1.0954451	9	6
Isacia_conceptionis	7	25.571429	2.5071327	30	22
Labrisomus_philippii	18	28.055556	6.9745101	38	15
Nexilosus_latifrons	4	23.000000	5.7154761	30	17
Paralabrax_humeralis	14	26.571429	8.4826908	44	15
Paralichthys_microps	4	32.500000	23.1012265	65	16
Paraxanthus_barbiger	12	6.250000	1.6583124	9	4
Patagonotothen_brevicauda	10	8.100000	1.3703203	10	6
Pinguipes_chilensis	81	28.370370	8.5197483	51	13
Prolatilus_jugularis	20	24.000000	6.9357958	39	10
Romaleon_setosum	17	8.058824	2.2768012	12	5
Scartichthys_spp.	176	13.670454	5.0927884	34	4
Schroederichthys_chilensis	3	51.666667	2.8867513	55	50
Sebastes_oculatus	12	12.333333	7.1137937	24	4
Semicossyphus_darwini	6	26.500000	6.4420494	36	20
Taliepus_sp.	27	9.444444	1.8045526	12	5

write.csv(Size_summary, "Size_summary2.csv")

Part 10: R Session Information

In this part, information about the operating system and R packages attached during the production of this document can be found

sessionInfo() %>% pander

R version 4.4.1 (2024-06-14 ucrt)

Platform: x86_64-w64-mingw32/x64

locale: LC_COLLATE=English_Australia.utf8, LC_CTYPE=English_Australia.utf8, LC_MONETARY=English_Australia.utf8, LC_NUMERIC=C and LC_TIME=English_Australia.utf8

attached base packages: stats, graphics, grDevices, utils, datasets, methods and base

other attached packages: mgcv(v.1.9-1), nlme(v.3.1-164), geosphere(v.1.5-18), MuMIn(v.1.48.4), MASS(v.7.3-60.2), corrplot(v.0.92), car(v.3.1-2), carData(v.3.0-5), knitr(v.1.47), performance(v.0.12.0), ggeffects(v.1.7.0), pander(v.0.6.5), patchwork(v.1.2.0), vegan(v.2.6-6.1), lattice(v.0.22-6), permute(v.0.9-7), summarytools(v.1.0.1), lubridate(v.1.9.3), forcats(v.1.0.0), stringr(v.1.5.1), dplyr(v.1.1.4), purrr(v.1.0.2), readr(v.2.1.5), tidyr(v.1.3.1), tibble(v.3.2.1), ggplot2(v.3.5.1) and tidyverse(v.2.0.0)

loaded via a namespace (and not attached): tidyselect(v.1.2.1), DHARMa(v.0.4.6), farver(v.2.1.2), fastmap(v.1.2.0), digest(v.0.6.35), timechange(v.0.3.0), lifecycle(v.1.0.4), cluster(v.2.1.6), magrittr(v.2.0.3), compiler(v.4.4.1), rlang(v.1.1.4), sass(v.0.4.9), tools(v.4.4.1), utf8(v.1.2.4), yaml(v.2.3.8), labeling(v.0.4.3), sp(v.2.1-4), plyr(v.1.8.9), abind(v.1.4-5), withr(v.3.0.0), grid(v.4.4.1), stats4(v.4.4.1), fansi(v.1.0.6), colorspace(v.2.1-0), scales(v.1.3.0), insight(v.0.20.1), cli(v.3.6.2), crayon(v.1.5.3), rmarkdown(v.2.27), ragg(v.1.3.2), generics(v.0.1.3), rstudioapi(v.0.16.0), reshape2(v.1.4.4), tzdb(v.0.4.0), minqa(v.1.2.7), cachem(v.1.1.0), splines(v.4.4.1), parallel(v.4.4.1), matrixStats(v.1.3.0), base64enc(v.0.1-3), vctrs(v.0.6.5), boot(v.1.3-30), Matrix(v.1.7-0), jsonlite(v.1.8.8), hms(v.1.1.3), rapportools(v.1.1), systemfonts(v.1.1.0), magick(v.2.8.3), jquerylib(v.0.1.4), glue(v.1.7.0), nloptr(v.2.1.0), codetools(v.0.2-20), stringi(v.1.8.4), gtable(v.0.3.5), lme4(v.1.1-35.4), munsell(v.0.5.1), pillar(v.1.9.0), htmltools(v.0.5.8.1), R6(v.2.5.1), textshaping(v.0.4.0), tcltk(v.4.4.1), evaluate(v.0.24.0), highr(v.0.11), backports(v.1.5.0), pryr(v.0.1.6), bslib(v.0.7.0), Rcpp(v.1.0.12), checkmate(v.2.3.1), xfun(v.0.44) and pkgconfig(v.2.0.3)

Environmental and local habitat variables as predictors of trophic interactions in subtidal rocky reefs along the SE Pacific coast

12-08-2024

Clear the environment

Load packages

Part 1: Importing and examining data

1a) Import data

1b) Inspect structure of data frames

Part 2: Manipulating datasets

2a) Environmental data

2b): Benthic data

2c): Fish length

2d) Species list (abundances and richness)

2e): Community composition

2f): Phlorotannins

Part 3: Latitudinal patterns in explanatory variables

3a): Explore latitudinal relationships

1) Temperature vs latitude

2) Salinity vs latitude

3) Oxygen vs latitude

4) Visibility vs latitude

5) PCOA vs latitude

6) Kelp vs latitude

7) Bare_rock vs latitude

8) Predatory fish abundance vs latitude

9) Herbivorous fish abundance vs latitude

10) Invertebrate predators abundance vs latitude

11) Species richness vs latitude

12) Fish length vs latitude

13) Phlorotannins vs latitude

Final plots

Part 4: Latitudinal patterns in consumption

4a): Predation vs latitude

1) Crabs 2h

2) Crabs 24h

3) Crabs 2 and 24 h prediction plot

4) Squidpops 2h

5) Squidpops 24h

3) Squidpops 2 and 24 h prediction plot

4b): Herbivory vs latitude

Part 5: Predation intensity: Effect of site-specific variables

5a): Standardisation and test autocorrelation

5b) Model selection: Crabs 2h

5c) Model selection: Crabs 24h

5e): Model selection: Squidpops 2h

5f) Model selection: Squidpops 24h

Part 6: Herbivory intensity: Effect of site specific variables

6a) Model selection

Part 7: Visualisation of best fitted predictors

7a): Crab consumption

7b): Squidpop consumption

7c): Seaweed consumption

Part 9: Supplementary tables

9a) Size estimations summary table

Part 10: R Session Information