# Load required libraries library(nlme) library(ggplot2) library(dplyr) library(scales) library(MuMIn) library(glmmTMB) library(corrplot) library(Hmisc) library(grid) library(gridGraphics) library(tidyr) library(stringr) library(vegan) library(ggplot2) library(ggnewscale) library(reshape) library(gridExtra) library(reshape2) library(tidyverse) library(ggpubr) #create colour palletes pal<-c("#76c68f", "grey", "#bd7ebe") pal2<-c("#76c68f", "#bd7ebe") pal3<- c("tomato4", "grey", "lightcoral") # Set working directories setwd("C:/Users/as21769/OneDrive - University of Essex/Documents/Essex/Phd/2023/New PPMR") PPMR_All_df <- read.csv("PPMR_Final_Df.csv") %>% unique() %>% as.data.frame() # Subset and rename relevant columns PPMR_All_df_sub <- PPMR_All_df %>% dplyr::select( ppmr, pred_weight_g, prey_ind_weight_g, prey_count, haul_id, stomach_id, year, Avg_SST, Mean_effort, Depth, meanSal, chl, Quarter, data, year_since, ICES, pred_taxa, prey_taxa, latitude ) %>% dplyr::rename( PPMR = ppmr, Pred_weight = pred_weight_g, Prey_weight = prey_ind_weight_g, Prey_count = prey_count, Year = year, Data_source = data, Pred_taxa = pred_taxa, Prey_taxa = prey_taxa, Latitude = latitude ) %>% unique() # DATA CLEANING # All PPMR data df_cleaned_PPMR <- PPMR_All_df_sub %>% filter(complete.cases(PPMR, Avg_SST), PPMR >= 1) %>% mutate( ICES = as.factor(ICES), Year = as.numeric(Year), Quarter = as.factor(Quarter), PPMR = as.numeric(PPMR), Log_PPMR = as.numeric(log10(PPMR)), Avg_SST = as.numeric(Avg_SST), Mean_effort = as.numeric(Mean_effort), year_since=(as.numeric(year_since)), Depth = as.numeric(Depth) ) %>% group_by(ICES) %>% filter(n() > 100) %>% ungroup() # subset of PPMR with FE df_cleaned_PPMR_FE <- PPMR_All_df_sub %>% filter(complete.cases(PPMR, Avg_SST, Mean_effort), PPMR >= 1) %>% mutate( ICES = as.factor(ICES), Year = as.numeric(Year), Quarter = as.factor(Quarter), PPMR = as.numeric(PPMR), Log_PPMR = as.numeric(log10(PPMR)), Avg_SST = as.numeric(Avg_SST), Mean_effort = as.numeric(Mean_effort), year_since=(as.numeric(year_since)), Depth = as.numeric(Depth) ) %>% group_by(ICES) %>% filter(n() > 100) %>% ungroup() # subset of All PPMR data with predator weight df_cleaned_PPMR_Pred <- df_cleaned_PPMR %>% filter(complete.cases(Pred_weight)) df_cleaned_PPMR_Pred <- df_cleaned_PPMR_Pred %>% # Aggregating data group_by(stomach_id) %>% dplyr::summarise( Year = Year, ICES = ICES, Pred_weight = mean(Pred_weight), Avg_SST = mean(Avg_SST), meanSal = mean(meanSal), chl = mean(chl), Data_source = Data_source, year_since = year_since, Quarter = Quarter, Pred_taxa = Pred_taxa, Prey_weight = mean(Prey_weight), Depth = mean(Depth) ) %>% unique() df_cleaned_PPMR_Pred <- df_cleaned_PPMR_Pred %>% mutate( log10_Pred_wg = as.numeric(log10(Pred_weight)), Avg_SST = as.numeric(Avg_SST), Prey_weight = as.numeric(Prey_weight), Log10_Prey_weight = as.numeric(log10(Prey_weight)) ) # subset of All PPMR data with prey weight df_cleaned_PPMR_Prey<- df_cleaned_PPMR %>% filter(complete.cases(Prey_weight)) # subset of PPMR data with FE and with predator weight df_cleaned_PPMR_Pred_FE <- df_cleaned_PPMR_FE %>% filter(complete.cases(Pred_weight)) df_cleaned_PPMR_Pred_FE <- df_cleaned_PPMR_Pred_FE %>% mutate( ICES = as.factor(ICES), Year = as.factor(Year), Quarter = as.factor(Quarter), Pred_weight = as.numeric(Pred_weight), Avg_SST = as.numeric(Avg_SST), log10_Pred_wg = log10(Pred_weight), Mean_effort = as.numeric(Mean_effort), year_since = as.numeric(year_since), Prey_weight = as.numeric(Prey_weight), Prey_weight_class = as.numeric(cut_number(Prey_weight, 3)), FE_class = as.numeric(cut_number(Mean_effort, 3)), Pred_taxa = as.factor(Pred_taxa) ) df_cleaned_PPMR_Pred_FE <- df_cleaned_PPMR_Pred_FE %>% # Group by stomach_id and calculate means group_by(stomach_id) %>% dplyr::summarize( ICES = first(ICES), log10_Pred_wg = mean(log10_Pred_wg), Avg_SST = mean(Avg_SST), Mean_effort = mean(Mean_effort), meanSal = mean(meanSal), chl = mean(chl), Data_source = Data_source, year_since = year_since, Quarter = Quarter, Prey_weight = mean(Prey_weight), FE_class = mean(FE_class), Pred_taxa = as.factor(Pred_taxa), Depth = mean(Depth) ) %>% distinct() # subset of PPMR data with FE and with prey weight df_cleaned_PPMR_Prey_FE<- df_cleaned_PPMR_FE %>% filter(complete.cases(Prey_weight)) # subset of PPMR data with prey weight and prey count df_cleaned_PPMR_Prey_Count<- df_cleaned_PPMR_Prey %>% filter(complete.cases(df_cleaned_PPMR$Prey_count)) # subset of PPMR data with FE, prey weight and prey count df_cleaned_PPMR_Prey_Count_FE<- df_cleaned_PPMR_Prey_FE %>% filter(complete.cases(Prey_count)) ##Correlation testing corr.variables<-as.data.frame(cbind(df_cleaned_PPMR$year_since, df_cleaned_PPMR$Avg_SST, df_cleaned_PPMR$chl, df_cleaned_PPMR$Mean_effort, df_cleaned_PPMR$meanSal, df_cleaned_PPMR$Depth, df_cleaned_PPMR$Quarter)) colnames(corr.variables)<-c("year_since", "Avg_SST", "chl", "Effort_mean", "meanSal", "Depth", "Quarter") cor<-cor(corr.variables,use="complete.obs") corrplot(cor) # PPMR vs SST PPMRLmmodel <- lme(Log_PPMR ~ Avg_SST, random = list(meanSal = ~1, chl = ~1 , ICES = ~ 1, year_since = ~1, Pred_taxa = ~1, Quarter = ~1 , stomach_id = ~1, Depth = ~1 ), data = df_cleaned_PPMR, na.action=na.exclude) PPMRLmmodel_poly <- lme(Log_PPMR ~ I(Avg_SST^2), random = list(meanSal = ~1, chl = ~1 , ICES = ~ 1, year_since = ~1, Quarter = ~1 , Pred_taxa = ~1 , stomach_id = ~1, Depth = ~1 ), data = df_cleaned_PPMR, na.action=na.exclude) #Model Selection AIC(PPMRLmmodel, PPMRLmmodel_poly) VarCorr(PPMRLmmodel) setwd("C:/Users/as21769/OneDrive - University of Essex/Documents/Essex/Phd/2023/PPMR chapter/Resubmission April 2025") saveRDS(PPMRLmmodel, "PPMRLmmodel.rdata") PPMRLmmodel<-readRDS("PPMRLmmodel.rdata") # Model Results summary(PPMRLmmodel) r.squaredGLMM(PPMRLmmodel) # Model Validation # Histogram df_cleaned_PPMR$residuals <- residuals(PPMRLmmodel) hist(df_cleaned_PPMR$residuals, breaks = 30, main = "") Histogram <- recordPlot() grid.echo(Histogram) Hist_plot <- grid.grab() # Residuals vs Fitted Plot plot(fitted(PPMRLmmodel), residuals(PPMRLmmodel), xlab = "Fitted Values", ylab = "Residuals") abline(h = 0, col = "red", lty = 2) residual_plot <- recordPlot() grid.echo(residual_plot) Resid_plot <- grid.grab() # Residuals Plot for each covariate covariates_cont <- c("Avg_SST", "meanSal", "chl", "year_since", "Depth") plot_list_cont <- list() plot_list_cont <- lapply(covariates_cont, function(covariate) { ggplot(df_cleaned_PPMR, aes_string(x = covariate, y = "residuals")) + geom_point(alpha = 0.5) + geom_hline(yintercept = 0, color = "red", linetype = "dashed") + labs(x = covariate, y = "Residuals") + theme_minimal() }) grob_list_cont <- lapply(plot_list_cont, ggplotGrob) covariates_Dis <- c("ICES", "Pred_taxa", "stomach_id", "Quarter") plot_list_Dis <- list() plot_list_Dis <- lapply(covariates_Dis, function(covariate) { ggplot(df_cleaned_PPMR, aes_string(x = covariate, y = "residuals")) + geom_boxplot(alpha = 0.5) + geom_hline(yintercept = 0, color = "red", linetype = "dashed") + labs(x = covariate, y = "Residuals") + theme_minimal() }) grob_list_Dis <- lapply(plot_list_Dis, ggplotGrob) all_plots <- list( Hist_plot, Resid_plot, grob_list_cont[[1]], grob_list_cont[[2]], grob_list_cont[[3]], grob_list_cont[[4]], grob_list_cont[[5]], grob_list_Dis[[1]], grob_list_Dis[[2]], grob_list_Dis[[3]], grob_list_Dis[[4]] ) labeled_plots <- mapply(function(plot, label) { arrangeGrob(plot, top = textGrob(label, x = unit(0, "npc"), y = unit(1, "npc"), just = c("left", "top"), gp = gpar(fontsize = 14, fontface = "bold"))) }, all_plots, letters[1:length(all_plots)], SIMPLIFY = FALSE) PPMR_diagnostics <- grid.arrange(grobs = labeled_plots, ncol = 4) # Predictions predictions_PPMR <- predict(PPMRLmmodel, newdata = df_cleaned_PPMR, level = 0, se = T) %>% as.data.frame() data_PPMR <- cbind(df_cleaned_PPMR, predictions_PPMR) # Plot PPMR vs SST p_SST<-ggplot() + geom_point(aes(x = df_cleaned_PPMR$Avg_SST, y = df_cleaned_PPMR$Log_PPMR), alpha = 0.03, size = 1)+ geom_smooth(aes(x = data_PPMR$Avg_SST, y = data_PPMR$fit), method = "lm", formula = y ~ x + I(x^2), color = "#015b58", linewidth = 2.5, se = TRUE) + labs(x= "Temperature (°C)", y="Log10 PPMR") + theme_classic() + theme(axis.text.x = element_text(size = 24, colour = "black"), axis.text.y = element_text(size = 24, colour = "black"), axis.title = element_text(face = "bold", size = 25)) + labs(subtitle = "b", size=10,fontface = "bold") #+ # stat_regline_equation( # data = data_PPMR, # aes(x = Avg_SST, y = fit, label = ..eq.label..), # label.y = 7.5, # label.x = 6 # ) p_SST model_data <- eval(PPMRLmmodel$call$data, environment(formula(PPMRLmmodel))) mf <- model.frame(formula(PPMRLmmodel), data = model_data, na.action = na.omit) mf<-unique(na.omit(mf)) plot_data <- ggplot_build(p_SST)$data[[2]] setwd("C:/Users/as21769/OneDrive - University of Essex/Documents/Essex/Phd/2023/PPMR chapter/FINAL PPMR PUBLICATION") write.csv(plot_data, "Figure 1 plot data.csv") fig1_df<-cbind(df_cleaned_PPMR$Avg_SST, df_cleaned_PPMR$Log_PPMR) write.csv(mf, "fig1_df.csv") length(df_cleaned_PPMR$Log_PPMR) #find the predicted change in PPMR SST_range<-range(df_cleaned_PPMR$Avg_SST) SST_data_PPMR <- as.data.frame(expand.grid(Avg_SST = SST_range)) predictions_PPMR_range <- predict(PPMRLmmodel, newdata = SST_data_PPMR, level = 0, se = T) %>% as.data.frame() per_dif<-((abs(predictions_PPMR_range[1,1]-predictions_PPMR_range[2,1]))/predictions_PPMR_range[1,1])*100 21.002144-3.930786 degreechange<-per_dif/17.07136 # PPMR vs SST and FE INTERACTION #original PPMRLmmodel_FE <- lme(Log_PPMR ~ Avg_SST+ Mean_effort + Avg_SST * Mean_effort, random = list(ICES = ~ 1, year_since = ~1, Quarter = ~1, chl = ~1, Pred_taxa = ~1, stomach_id = ~1, meanSal = ~1, Depth = ~1 ), data = df_cleaned_PPMR_FE, na.action=na.exclude) PPMRLmmodel_FE_poly <- lme(Log_PPMR ~ I(Avg_SST^2)+ I(Mean_effort^2) + I(Avg_SST^2) * I(Mean_effort^2), random = list( ICES = ~ 1, year_since = ~1, Quarter = ~1, chl = ~1, Pred_taxa = ~1, stomach_id = ~1, meanSal = ~1, Depth = ~1), data = df_cleaned_PPMR_FE, na.action=na.exclude) AIC(PPMRLmmodel_FE, PPMRLmmodel_FE_poly) unique(df_cleaned_PPMR_FE) df_clean <- df %>% filter(!is.na(x)) setwd("C:/Users/as21769/OneDrive - University of Essex/Documents/Essex/Phd/2023/PPMR chapter/Resubmission April 2025") saveRDS(PPMRLmmodel_FE, "PPMRLmmodel_FE.rdata") PPMRLmmodel_FE<-readRDS("PPMRLmmodel_FE.rdata") # Model Results summary(PPMRLmmodel_FE) r.squaredGLMM(PPMRLmmodel_FE) # Model Validation # Histogram residuals <- resid(PPMRLmmodel_FE) hist(residuals, breaks = 30, main = "") Histogram <- recordPlot() grid.echo(Histogram) Hist_plot <- grid.grab() # Fitted vs Residual Plot plot(fitted(PPMRLmmodel_FE), residuals(PPMRLmmodel_FE), xlab = "Fitted Values", ylab = "Residuals") abline(h = 0, col = "red", lty = 2) residual_plot <- recordPlot() grid.echo(residual_plot) Resid_plot <- grid.grab() # Residuals Plot for each covariate df_cleaned_PPMR_FE$residuals <- residuals(PPMRLmmodel_FE) covariates_cont <- c("Avg_SST", "Mean_effort", "meanSal", "chl", "year_since","Depth") plot_list_cont <- list() plot_list_cont <- lapply(covariates_cont, function(covariate) { ggplot(df_cleaned_PPMR_FE, aes_string(x = covariate, y = "residuals")) + geom_point(alpha = 0.5) + geom_hline(yintercept = 0, color = "red", linetype = "dashed") + labs(x = covariate, y = "Residuals") + theme_minimal() }) grob_list_cont <- lapply(plot_list_cont, ggplotGrob) covariates_Dis <- c("ICES", "Pred_taxa", "stomach_id", "Quarter") plot_list_Dis <- list() plot_list_Dis <- lapply(covariates_Dis, function(covariate) { ggplot(df_cleaned_PPMR_FE, aes_string(x = covariate, y = "residuals")) + geom_boxplot(alpha = 0.5) + geom_hline(yintercept = 0, color = "red", linetype = "dashed") + labs(x = covariate, y = "Residuals") + theme_minimal() }) grob_list_Dis <- lapply(plot_list_Dis, ggplotGrob) all_plots_FE <- list( Hist_plot, Resid_plot, grob_list_cont[[1]], grob_list_cont[[2]], grob_list_cont[[3]], grob_list_cont[[4]], grob_list_cont[[5]], grob_list_cont[[6]], grob_list_Dis[[1]], grob_list_Dis[[2]], grob_list_Dis[[3]], grob_list_Dis[[4]] ) labeled_plots_FE <- mapply(function(plot, label) { arrangeGrob(plot, top = textGrob(label, x = unit(0, "npc"), y = unit(1, "npc"), just = c("left", "top"), gp = gpar(fontsize = 14, fontface = "bold"))) }, all_plots_FE, letters[1:length(all_plots_FE)], SIMPLIFY = FALSE) PPMR_diagnostics_FE <- grid.arrange(grobs = labeled_plots_FE, ncol = 4) #Find the qartiles of SST for predictions SST_quantiles<-as.numeric(quantile(df_cleaned_PPMR_FE$Avg_SST)) # Binning Fishing Effort df_cleaned_PPMR_FE <- df_cleaned_PPMR_FE %>% mutate(FEE = as.numeric(cut_number(Mean_effort, 3))) # Calculate Medians for predictions FE_medians <- df_cleaned_PPMR_FE %>% group_by(FEE) %>% dplyr::summarize(Median_Effort = median(Mean_effort, na.rm = TRUE)) %>% pull(Median_Effort) FE_medians<-na.omit(FE_medians) new_data_PPMR <- as.data.frame(expand.grid(Avg_SST = SST_quantiles, Mean_effort = FE_medians)) new_data_PPMR$Avg_SSTx2<-(new_data_PPMR$Avg_SST)^2 new_data_PPMR$Mean_effortx2<-(new_data_PPMR$Mean_effort)^2 predictions_PPMR <- predict(PPMRLmmodel_FE, newdata = new_data_PPMR, level = 0, se = TRUE) %>% as.data.frame() dataInt_PPMR <- cbind(new_data_PPMR, predictions_PPMR) dataInt_PPMR<- dataInt_PPMR %>% mutate(FEE = case_when(Mean_effort == FE_medians[1] ~ "1", Mean_effort == FE_medians[2] ~ "2", Mean_effort == FE_medians[3] ~ "3" )) df_cleaned_PPMR_FE$FEE<-as.factor(df_cleaned_PPMR_FE$FEE) dataInt_PPMR$FEE<-as.factor(dataInt_PPMR$FEE) dataInt_PPMR$Mean_effort<-as.numeric(dataInt_PPMR$Mean_effort) dataInt_PPMR$Avg_SST<-as.numeric(dataInt_PPMR$Avg_SST) dataInt_PPMR$fit<-as.numeric(dataInt_PPMR$fit) df_cleaned_PPMR_FE$Avg_SST<-as.numeric(df_cleaned_PPMR_FE$Avg_SST) df_cleaned_PPMR_FE$Log_PPMR<-as.numeric(df_cleaned_PPMR_FE$Log_PPMR) df_cleaned_PPMR_FE$Mean_effort<-as.numeric(df_cleaned_PPMR_FE$Mean_effort) p_FE<-ggplot() + geom_point(data = df_cleaned_PPMR_FE, mapping = aes(x=Avg_SST, y=Log_PPMR, color=FEE), group=df_cleaned_PPMR_FE$FEE, alpha = 0.025, size = 1) + geom_ribbon(data= dataInt_PPMR, aes(x = Avg_SST, ymin= fit - se.fit, ymax = fit +se.fit, color = FEE), fill = "grey", alpha = 0.25) + geom_line(data= dataInt_PPMR, aes(x = Avg_SST, y = fit, color = FEE), linewidth = 1.5) + scale_color_manual(labels=c('Low', 'Medium', 'High'), values = pal) + scale_fill_manual(labels=c('Low', 'Medium', 'High'), values = pal) + labs(x= "Temperature (°C)", y="Log10 PPMR", color="Mean Fishing\nEffort\n(hours per year)", fill="Mean Fishing\nEffort\n(hours per year)") + theme_classic() + theme(axis.text.x=element_text(size=16,colour="black")) + theme(axis.text.y=element_text(size=16,colour="black")) + theme(axis.title=element_text(face="bold", size=18)) + theme(legend.position.inside = c(0.8,0.9), legend.text = element_text(size=12),legend.title = element_text(size=14))+ labs(subtitle = "d", size=10,fontface = "bold") p_FE model_data <- eval(PPMRLmmodel_FE$call$data, environment(formula(PPMRLmmodel_FE))) mf <- model.frame(formula(PPMRLmmodel_FE), data = model_data, na.action = na.omit) mf<-unique(na.omit(mf)) plot_data <- ggplot_build(p_FE)$data[[2]] setwd("C:/Users/as21769/OneDrive - University of Essex/Documents/Essex/Phd/2023/PPMR chapter/FINAL PPMR PUBLICATION") write.csv(plot_data, "Figure 1 plot data_FE.csv") write.csv(mf, "fig1_df_FE.csv") p_FE<-ggplot() + geom_smooth(data= dataInt_PPMR, aes(x = Avg_SST, y= fit, color = FEE), fill = "grey", alpha = 0.25, method = "lm", size = 2) + geom_point(data = df_cleaned_PPMR_FE, mapping = aes(x=Avg_SST, y=Log_PPMR, color=FEE), group=df_cleaned_PPMR_FE$FEE, alpha = 0.025) + scale_color_manual(labels=c('Low', 'Medium', 'High'), values = pal) + scale_fill_manual(labels=c('Low', 'Medium', 'High'), values = pal) + labs(x= "Temperature (°C)", y="Log10 PPMR", color="Mean Fishing\nEffort\n(hours per year)", fill="Mean Fishing\nEffort\n(hours per year)") + theme_classic() + theme(axis.text.x=element_text(size=16,colour="black")) + theme(axis.text.y=element_text(size=16,colour="black")) + theme(axis.title=element_text(face="bold", size=18)) + theme(legend.position.inside = c(0.8,0.9), legend.text = element_text(size=12),legend.title = element_text(size=14))+ stat_regline_equation( data = dataInt_PPMR, aes(x = Avg_SST, y = fit, label = ..eq.label.., color = FEE), label.y = c(5.1, 5.6, 6), label.x = 6 ) p_FE #find the predicted change in PPMR SST_FE_data_PPMR <- as.data.frame(expand.grid(Avg_SST = SST_range, Mean_effort = FE_medians)) predictions_PPMR_FE_range <- predict(PPMRLmmodel_FE, newdata = SST_FE_data_PPMR, level = 0, se = T) %>% as.data.frame() predictions_PPMR_FE_range_df <- cbind(SST_FE_data_PPMR, predictions_PPMR_FE_range) per_dif_FE_low<-((abs(predictions_PPMR_FE_range_df[1,3]-predictions_PPMR_FE_range_df[2,3]))/predictions_PPMR_FE_range_df[1,3])*100 per_dif_FE_high<-((abs(predictions_PPMR_FE_range_df[5,3]-predictions_PPMR_FE_range_df[6,3]))/predictions_PPMR_FE_range_df[5,3])*100 degreechange_FE_low<-per_dif_FE_low/17.07136 degreechange_FE_high<-per_dif_FE_high/17.07136 ##PREDATOR BODY MASS Pred_wgLmmodel <- lme(log10_Pred_wg ~ Avg_SST, random = list(ICES = ~ 1, year_since = ~1, Quarter = ~1, chl = ~1, Pred_taxa = ~1, stomach_id = ~1, meanSal = ~1, Depth = ~1), data = df_cleaned_PPMR_Pred, na.action = na.exclude) Pred_wgLmmodel_poly <- lme(log10_Pred_wg ~ I(Avg_SST^2), random = list(ICES = ~ 1, year_since = ~1, Quarter = ~1, chl = ~1, Pred_taxa = ~1, stomach_id = ~1, meanSal = ~1, Depth = ~1), data = df_cleaned_PPMR_Pred, na.action = na.exclude) AIC(Pred_wgLmmodel, Pred_wgLmmodel_poly) #Model Results setwd("C:/Users/as21769/OneDrive - University of Essex/Documents/Essex/Phd/2023/PPMR chapter/Resubmission April 2025") saveRDS(Pred_wgLmmodel, "Pred_wgLmmodel.rdata") Pred_wgLmmodel<-readRDS("Pred_wgLmmodel.rdata") summary(Pred_wgLmmodel) r.squaredGLMM(Pred_wgLmmodel) # Predictions and plotting predictions_pred_wg <- predict(Pred_wgLmmodel, newdata = df_cleaned_PPMR_Pred, level = 0, se.fit = TRUE, type = "link") %>% as.data.frame() data_pred_wg <- cbind(df_cleaned_PPMR_Pred, predictions_pred_wg) data_pred_wg$fit <- as.numeric(data_pred_wg$fit) data_pred_wg$Avg_SST <- as.numeric(data_pred_wg$Avg_SST) y.expression.pred <- expression(bold(Log[10] ~ Predator ~ Body ~ Mass)) p_Pred_wg_all <- ggplot() + geom_point(data = df_cleaned_PPMR_Pred, aes(x = Avg_SST, y = log10_Pred_wg), alpha = 0.01, size = 1) + geom_smooth(data = data_pred_wg, aes(x = Avg_SST, y = fit), method = "lm", linewidth = 1.5, color = "indianred3", fullrange = T) + labs(x = "Temperature (°C)", y = y.expression.pred) + theme_classic() + theme( axis.text.x = element_text(size = 24, colour = "black"), axis.text.y = element_text(size = 24, colour = "black"), axis.title = element_text(face = "bold", size = 25), legend.position = "none" ) + labs(subtitle = "a", size = 10, fontface = "bold") + theme(plot.subtitle = element_text(size = 26, face = "bold", family = "arial"))#+ # stat_regline_equation( # data = data_pred_wg, # aes(x = Avg_SST, y = fit, label = ..eq.label..), # label.y = c(5.1, 5.6, 6), # label.x = 6 # ) p_Pred_wg_all setwd("C:/Users/as21769/OneDrive - University of Essex/Documents/Essex/Phd/2023/PPMR chapter/FINAL PPMR PUBLICATION") model_data <- eval(Pred_wgLmmodel$call$data, environment(formula(Pred_wgLmmodel))) mf <- model.frame(formula(Pred_wgLmmodel), data = model_data, na.action = na.omit) mf<-unique(na.omit(mf)) write.csv(mf, "Pred_wgLmmodel.csv") plot_data <- ggplot_build(p_Pred_wg_all)$data[[2]] write.csv(plot_data, "Figure 1 Pred_wgLmmodel.csv") ## Predator Interaction Pred_wgLmmodel_FE <- lme(log10_Pred_wg ~ Avg_SST + Mean_effort + Avg_SST:Mean_effort, random = list(ICES = ~ 1, year_since = ~1, Quarter = ~1, chl = ~1, meanSal = ~1, Depth = ~1, Pred_taxa = ~1 ), data = df_cleaned_PPMR_Pred_FE, na.action = na.exclude) Pred_wgLmmodel_poly_FE <- lme(I(log10_Pred_wg^2) ~ I(Avg_SST^2) + I(Mean_effort^2) + I(Avg_SST^2):I(Mean_effort^2), random = list(ICES = ~ 1, year_since = ~1, Quarter = ~1, chl = ~1, meanSal = ~1, Depth = ~1 , Pred_taxa = ~1 ), data = df_cleaned_PPMR_Pred_FE, na.action = na.exclude) AIC(Pred_wgLmmodel_FE, Pred_wgLmmodel_poly_FE) setwd("C:/Users/as21769/OneDrive - University of Essex/Documents/Essex/Phd/2023/PPMR chapter/Resubmission April 2025") saveRDS(Pred_wgLmmodel_FE, "Pred_wgLmmodel_FE.rdata") Pred_wgLmmodel_FE<-readRDS("Pred_wgLmmodel_FE.rdata") # Model Results summary(Pred_wgLmmodel_FE) r.squaredGLMM(Pred_wgLmmodel_FE) # Predictions new_data_Pred <- expand.grid(Avg_SST = SST_quantiles, Mean_effort = FE_medians) predictions_PredInt <- predict(Pred_wgLmmodel_FE, newdata = new_data_Pred, level = 0, se = TRUE) %>% as.data.frame() dataInt_Pred <- cbind(new_data_Pred, predictions_PredInt) dataInt_Pred<- dataInt_Pred %>% mutate(FE_class = case_when(Mean_effort == FE_medians[1] ~ "1", Mean_effort == FE_medians[2] ~ "2", Mean_effort == FE_medians[3] ~ "3" )) # Plot p_Pred_wg_Int <- ggplot() + geom_point(data = df_cleaned_PPMR_Pred_FE, aes(x = Avg_SST, y = log10_Pred_wg, color = as.factor(FE_class)), alpha = 0.25, size = 1) + geom_smooth(method = "lm", aes(x = Avg_SST, y = fit, color = as.factor(FE_class)), data = subset(dataInt_Pred, FE_class %in% c("1", "3")), size = 1.5, se = TRUE) + scale_color_manual(labels = c('Low', 'Medium', 'High'), values = c("#76c68f", "grey", "#bd7ebe")) + labs( x = "Temperature (°C)", y = expression(bold(Log[10] ~ Predator ~ Body ~ Mass)), color = "Mean Fishing Effort" ) + theme_classic() + theme( axis.text.x = element_text(size = 24, colour = "black"), axis.text.y = element_text(size = 24, colour = "black"), axis.title = element_text(face = "bold", size = 25), legend.position = "none" , plot.subtitle = element_text(size = 26, face = "bold", family = "arial") )+ labs(subtitle = "a", size = 10, fontface = "bold")#+ # stat_regline_equation( # data = dataInt_Pred, # aes(x = Avg_SST, y = fit, label = ..eq.label.., color = FE_class), # label.y = c(5.1, 5.6, 6), # label.x = 6 # ) p_Pred_wg_Int setwd("C:/Users/as21769/OneDrive - University of Essex/Documents/Essex/Phd/2023/PPMR chapter/FINAL PPMR PUBLICATION") model_data <- eval(Pred_wgLmmodel_FE$call$data, environment(formula(Pred_wgLmmodel_FE))) mf <- model.frame(formula(Pred_wgLmmodel_FE), data = model_data, na.action = na.omit) mf<-unique(na.omit(mf)) write.csv(mf, "Pred_wgLmmodel_FE.csv") plot_data <- ggplot_build(p_Pred_wg_Int)$data[[2]] write.csv(plot_data, "Figure 1 Pred_wgLmmodel_FE.csv") ## Prey Body mass all df_cleaned_PPMR_Prey$Prey_weight_class<-as.numeric(cut_number(df_cleaned_PPMR_Prey$Prey_weight,3)) one<-subset(df_cleaned_PPMR_Prey, Prey_weight_class == 1) one<-one[complete.cases(one$Prey_weight), ] median_one<-mean(one$Prey_weight) two<-subset(df_cleaned_PPMR_Prey, Prey_weight_class == 2) two<-two[complete.cases(two$Prey_weight), ] median_two<-median(two$Prey_weight) three<-subset(df_cleaned_PPMR_Prey, Prey_weight_class == 3) three<-three[complete.cases(three$Prey_weight), ] median_three<-median(three$Prey_weight) df_cleaned_PPMR_Prey$log10_Prey_wg<-log10(as.numeric(df_cleaned_PPMR_Prey$Prey_weight)) df_cleaned_PPMR_Prey$Avg_SST<-as.numeric(df_cleaned_PPMR_Prey$Avg_SST) df_cleaned_PPMR_Prey$meanSal<-as.numeric(df_cleaned_PPMR_Prey$meanSal) df_cleaned_PPMR_Prey$chl<-as.numeric(df_cleaned_PPMR_Prey$chl) df_cleaned_PPMR_Prey$year_since<-as.numeric(df_cleaned_PPMR_Prey$year_since) df_cleaned_PPMR_Prey$Prey_weight_class<-as.numeric(df_cleaned_PPMR_Prey$Prey_weight_class) df_cleaned_PPMR_Prey[sapply(df_cleaned_PPMR_Prey, is.infinite)] <- NA Prey_wgLmmodel <- lme(log10_Prey_wg ~ Avg_SST, random = list(ICES = ~ 1, year_since = ~1, Quarter = ~1, chl = ~1, meanSal = ~1, Depth = ~1, Prey_weight_class= ~1 ), data = df_cleaned_PPMR_Prey, na.action = na.exclude) Prey_wgLmmodel_poly <- lme(log10_Prey_wg ~ I(Avg_SST^2), random = list(ICES = ~ 1, year_since = ~1, Quarter = ~1, chl = ~1, meanSal = ~1, Depth = ~1, Prey_weight_class= ~1 ), data = df_cleaned_PPMR_Prey, na.action = na.exclude) AIC(Prey_wgLmmodel, Prey_wgLmmodel_poly) setwd("C:/Users/as21769/OneDrive - University of Essex/Documents/Essex/Phd/2023/PPMR chapter/Resubmission April 2025") saveRDS(Prey_wgLmmodel_poly, "Prey_wgLmmodel_poly.rdata") Prey_wgLmmodel_poly<-readRDS("Prey_wgLmmodel_poly.rdata") #Model Results summary(Prey_wgLmmodel_poly) r.squaredGLMM(Prey_wgLmmodel_poly) # Predictions and plotting predictions_Prey_wg <- predict(Prey_wgLmmodel_poly, newdata = df_cleaned_PPMR_Prey, level = 0, se.fit = TRUE, type = "link") %>% as.data.frame() data_Prey_wg <- cbind(df_cleaned_PPMR_Prey, predictions_Prey_wg) data_Prey_wg$fit <- as.numeric(data_Prey_wg$fit) data_Prey_wg$Avg_SST <- as.numeric(data_Prey_wg$Avg_SST) y.expression.Prey <- expression(bold(Log[10] ~ Prey ~ Body ~ Mass)) p_Prey_wg_all <- ggplot() + geom_point(data = df_cleaned_PPMR_Prey, aes(x = Avg_SST, y = log10_Prey_wg), alpha = 0.01, size = 1) + geom_smooth(data = data_Prey_wg, aes(x = Avg_SST, y = fit), method = "lm", linewidth = 1.5, formula = y ~ I(x^2), se = T, color = "indianred3", fullrange = T) + labs(x = "Temperature (°C)", y = y.expression.Prey) + theme_classic() + theme( axis.text.x = element_text(size = 24, colour = "black"), axis.text.y = element_text(size = 24, colour = "black"), axis.title = element_text(face = "bold", size = 25), legend.position = "none" ) + labs(subtitle = "b", size = 10, fontface = "bold") + theme(plot.subtitle = element_text(size = 26, face = "bold", family = "arial")) #+ #stat_regline_equation( #data = data_Prey_wg, #aes(x = Avg_SST, y = fit, label = ..eq.label..), #label.y = 3, #label.x = 6, #formula = y ~ poly(x, 2, raw = TRUE) #) p_Prey_wg_all setwd("C:/Users/as21769/OneDrive - University of Essex/Documents/Essex/Phd/2023/PPMR chapter/FINAL PPMR PUBLICATION") model_data <- eval(Prey_wgLmmodel_poly$call$data, environment(formula(Prey_wgLmmodel_poly))) mf <- model.frame(formula(Prey_wgLmmodel_poly), data = model_data, na.action = na.omit) mf<-unique(na.omit(mf)) write.csv(mf, "Prey_wgLmmodel_poly.csv") plot_data <- ggplot_build(p_Prey_wg_all)$data[[2]] write.csv(plot_data, "Figure 1 Prey_wgLmmodel_poly.csv") ##Prey body mass interaction df_cleaned_PPMR_Prey_FE$Prey_weight_class<-as.numeric(cut_number(df_cleaned_PPMR_Prey_FE$Prey_weight,3)) df_cleaned_PPMR_Prey_FE$log10_Prey_wg<-log10(as.numeric(df_cleaned_PPMR_Prey_FE$Prey_weight)) df_cleaned_PPMR_Prey_FE$Avg_SST<-as.numeric(df_cleaned_PPMR_Prey_FE$Avg_SST) df_cleaned_PPMR_Prey_FE$meanSal<-as.numeric(df_cleaned_PPMR_Prey_FE$meanSal) df_cleaned_PPMR_Prey_FE$chl<-as.numeric(df_cleaned_PPMR_Prey_FE$chl) df_cleaned_PPMR_Prey_FE$year_since<-as.numeric(df_cleaned_PPMR_Prey_FE$year_since) df_cleaned_PPMR_Prey_FE$Prey_weight_class<-as.numeric(df_cleaned_PPMR_Prey_FE$Prey_weight_class) df_cleaned_PPMR_Prey_FE[sapply(df_cleaned_PPMR_Prey_FE, is.infinite)] <- NA df_cleaned_PPMR_Prey_FE$FE_class = as.numeric(cut_number(df_cleaned_PPMR_Prey_FE$Mean_effort, 3)) df_cleaned_PPMR_Prey_FE$log10_Prey_wg<-as.numeric(df_cleaned_PPMR_Prey_FE$log10_Prey_wg) df_cleaned_PPMR_Prey_FE$Avg_SST<-as.numeric(df_cleaned_PPMR_Prey_FE$Avg_SST) df_cleaned_PPMR_Prey_FE$Mean_effort<-as.numeric(df_cleaned_PPMR_Prey_FE$Mean_effort) Prey_wgLmmodel_FE <- lme(log10_Prey_wg ~ Avg_SST+ Mean_effort + Avg_SST*Mean_effort, random = list(ICES = ~ 1, year_since = ~1, Quarter = ~1, chl = ~1, meanSal = ~1, Depth = ~1, Prey_weight_class= ~1 ), data = df_cleaned_PPMR_Prey_FE, na.action = na.exclude) Prey_wgLmmodel_poly_FE <- lme(log10_Prey_wg ~ I(Avg_SST^2)+ I(Mean_effort^2) + I(Avg_SST^2)*I(Mean_effort^2), random = list(ICES = ~ 1, year_since = ~1, Quarter = ~1, chl = ~1, meanSal = ~1, Depth = ~1, Prey_weight_class= ~1 ), data = df_cleaned_PPMR_Prey_FE, na.action = na.exclude) AIC(Prey_wgLmmodel_FE, Prey_wgLmmodel_poly_FE) setwd("C:/Users/as21769/OneDrive - University of Essex/Documents/Essex/Phd/2023/PPMR chapter/Resubmission April 2025") saveRDS(Prey_wgLmmodel_FE, "Prey_wgLmmodel_FE.rdata") Prey_wgLmmodel_FE<-readRDS("Prey_wgLmmodel_FE.rdata") #Model Results summary(Prey_wgLmmodel_FE) r.squaredGLMM(Prey_wgLmmodel_FE) #Predict new_data_Prey_wg <- expand.grid(Avg_SST = SST_quantiles, Mean_effort = FE_medians) predictions_Prey_wgInt<-predict(Prey_wgLmmodel_FE, newdata = new_data_Prey_wg, level = 0, se.fit = TRUE, type = "link") %>% as.data.frame() dataInt_Prey_wg<-as.data.frame(cbind(new_data_Prey_wg, predictions_Prey_wgInt)) dataInt_Prey_wg$Mean_effort<-as.factor(dataInt_Prey_wg$Mean_effort) df_cleaned_PPMR_Prey_FE$Mean_effort<-as.numeric(df_cleaned_PPMR_Prey_FE$Mean_effort) dataInt_Prey_wg$Avg_SST<-as.numeric(dataInt_Prey_wg$Avg_SST) dataInt_Prey_wg$fit<-as.numeric(dataInt_Prey_wg$fit) dataInt_Prey_wg<- dataInt_Prey_wg %>% mutate(FE_class = case_when(Mean_effort == FE_medians[1] ~ "1", Mean_effort == FE_medians[2] ~ "2", Mean_effort== FE_medians[3] ~ "3" )) df_cleaned_PPMR_Prey_FE$FE_class<-as.factor(df_cleaned_PPMR_Prey_FE$FE_class) dataInt_Prey_wg$Avg_SST<-as.numeric(dataInt_Prey_wg$Avg_SST) dataInt_Prey_wg$Mean_effort<-as.factor(dataInt_Prey_wg$Mean_effort) df_cleaned_PPMR_Prey_FE$Avg_SST<-as.numeric(df_cleaned_PPMR_Prey_FE$Avg_SST) p_Prey_wg_Int <- ggplot() + geom_point(data = df_cleaned_PPMR_Prey_FE, aes(x = Avg_SST, y = log10_Prey_wg, color = as.factor(FE_class)), alpha = 0.25, size = 1) + geom_smooth(method = "lm", aes(x = Avg_SST, y = fit, color = as.factor(FE_class)), data = subset(dataInt_Prey_wg, FE_class %in% c("1", "3")), size = 1.5, se = TRUE) + scale_color_manual(labels = c('Low', 'Medium', 'High'), values = c("#76c68f", "grey", "#bd7ebe")) + labs( x = "Temperature (°C)", y = expression(bold(Log[10] ~ Prey ~ Body ~ Mass)), color = "Mean Fishing Effort" ) + theme_classic() + theme( axis.text.x = element_text(size = 24, colour = "black"), axis.text.y = element_text(size = 24, colour = "black"), axis.title = element_text(face = "bold", size = 25), plot.subtitle = element_text(size = 26, face = "bold", family = "arial") )+ theme(legend.position = "none")+ labs(subtitle = "b", size = 10, fontface = "bold")#+ # stat_regline_equation( # data = dataInt_Prey_wg, # aes(x = Avg_SST, y = fit, label = ..eq.label.., color = FE_class), # label.y = c(2, 2.5, 3), # label.x = 6 # ) p_Prey_wg_Int setwd("C:/Users/as21769/OneDrive - University of Essex/Documents/Essex/Phd/2023/PPMR chapter/FINAL PPMR PUBLICATION") model_data <- eval(Prey_wgLmmodel_FE$call$data, environment(formula(Prey_wgLmmodel_FE))) mf <- model.frame(formula(Prey_wgLmmodel_FE), data = model_data, na.action = na.omit) mf<-unique(na.omit(mf)) write.csv(mf, "Prey_wgLmmodel_FE.csv") plot_data <- ggplot_build(p_Prey_wg_Int)$data[[2]] write.csv(plot_data, "Figure 1 Prey_wgLmmodel_FE.csv") ##Prey Count Prey_count_df<- df_cleaned_PPMR %>% group_by(stomach_id) %>% dplyr::summarize(Prey_count = sum(Prey_count), Avg_SST = Avg_SST, Mean_Prey_weight = mean(Prey_weight), Mean_effort= mean(Mean_effort), Pred_taxa = Pred_taxa, chl = mean(chl), year_since = year_since, ICES = ICES, Quarter = Quarter, meanSal = mean(meanSal), Depth = mean(Depth)) Prey_count_df$Mean_Prey_weight_class<-as.factor(cut_number(Prey_count_df$Mean_Prey_weight,3)) Prey_count_df$log10_Prey_count<-as.numeric(log10(Prey_count_df$Prey_count)) Prey_count_df<-as.data.frame(Prey_count_df) Prey_count_df$Avg_SST<-as.numeric(Prey_count_df$Avg_SST) Prey_count_df<-na.omit(Prey_count_df) Prey_count_df$FE_class<-as.numeric(cut_number(Prey_count_df$Mean_effort,3)) Prey_count_df<- Prey_count_df %>% mutate(Prey_weight_class = case_when(Mean_Prey_weight_class == "[1.85e-05,0.217]" ~ "1", Mean_Prey_weight_class == "(0.217,1.5]" ~ "2", Mean_Prey_weight_class== "(1.5,2.14e+03]" ~ "3" )) Prey_count_df$Prey_weight_class<-as.numeric(Prey_count_df$Prey_weight_class) Prey_countLmmodel <- lme(log10_Prey_count ~ Avg_SST+ Avg_SST*Prey_weight_class, random = list(ICES = ~ 1, year_since = ~1, Quarter = ~1, chl = ~1, meanSal = ~1, Depth = ~1, Pred_taxa = ~1 ), data = Prey_count_df, na.action = na.exclude) Prey_countLmmodel_poly <- lme(I(log10_Prey_count^2) ~ I(Avg_SST^2)+ I(Avg_SST^2)*I(Prey_weight_class^2), random = list(ICES = ~ 1, year_since = ~1, Quarter = ~1, chl = ~1, meanSal = ~1, Depth = ~1, Pred_taxa = ~1 ), data = Prey_count_df, na.action = na.exclude) AIC(Prey_countLmmodel, Prey_countLmmodel_poly) setwd("C:/Users/as21769/OneDrive - University of Essex/Documents/Essex/Phd/2023/PPMR chapter/Resubmission April 2025") saveRDS(Prey_countLmmodel, "Prey_countLmmodel.rdata") Prey_countLmmodel<-readRDS("Prey_countLmmodel.rdata") unique(df_cleaned_PPMR$stomach_id) #Model Results summary(Prey_countLmmodel) r.squaredGLMM(Prey_countLmmodel) #Model Predictions predictions_count<-predict(Prey_countLmmodel, newdata = Prey_count_df, level = 0, se.fit = TRUE, type = "link") %>% as.data.frame() data_count<-cbind(Prey_count_df, predictions_count) data_count$fit<-as.numeric(data_count$fit) data_count$Avg_SST<-as.numeric(data_count$Avg_SST) Prey_count_df$Mean_Prey_weight_class<-as.factor(Prey_count_df$Mean_Prey_weight_class) data_count$Mean_Prey_weight_class<-as.factor(data_count$Mean_Prey_weight_class) y.expression.count <- expression(bold(Log[10] ~ Mean ~ Prey ~ Count)) Prey_count_plot<-ggplot() + geom_point(aes(y=Prey_count_df$log10_Prey_count, x= Prey_count_df$Avg_SST,colour=Prey_count_df$Mean_Prey_weight_class),group=Prey_count_df$Mean_Prey_weight_class, alpha = 0.35, size = 1) + geom_smooth(method = "lm", aes(x = Avg_SST, y = fit, color = as.factor(Mean_Prey_weight_class)), data = data_count, size = 1.5, se = TRUE) + scale_color_manual(labels=c('Small', 'Medium', 'Large'), values = pal3) + scale_fill_manual(values = pal3) + labs(x= "Temperature (°C)", y=y.expression.count, color="Prey\nweight class", fill="Prey weight\nclass") + theme_classic() + theme(axis.text.x=element_text(size=24,colour="black")) + theme(axis.text.y=element_text(size=24,colour="black")) + theme(axis.title=element_text(face="bold", size=25)) + theme(legend.title =element_text(face="bold", size=18)) + theme(legend.text =element_text(size=18))+ theme(legend.position = "none")+ labs(subtitle = "c", size=10,fontface = "bold") + theme(plot.subtitle = element_text(size = 26, face = "bold",family="arial"))#+ # stat_regline_equation( # data = data_count, # aes(x = Avg_SST, y = fit, label = ..eq.label.., color =Mean_Prey_weight_class), # label.y = c(3.9, 4.2, 4.5), # label.x = 12.2 # ) Prey_count_plot setwd("C:/Users/as21769/OneDrive - University of Essex/Documents/Essex/Phd/2023/PPMR chapter/FINAL PPMR PUBLICATION") model_data <- eval(Prey_countLmmodel$call$data, environment(formula(Prey_countLmmodel))) mf <- model.frame(formula(Prey_countLmmodel), data = model_data, na.action = na.omit) mf<-unique(na.omit(mf)) write.csv(mf, "Prey_countLmmodel.csv") plot_data <- ggplot_build(Prey_count_plot)$data[[2]] write.csv(plot_data, "Figure 1 Prey_countLmmodel.csv") #Prey Count Interaction Prey_count_df_FE <- Prey_count_df %>% filter(complete.cases(Mean_effort)) Prey_count_df_FE<- Prey_count_df_FE %>% mutate(Prey_weight_class = case_when(Mean_Prey_weight_class == "[1.85e-05,0.217]" ~ "1", Mean_Prey_weight_class == "(0.217,1.5]" ~ "2", Mean_Prey_weight_class== "(1.5,2.14e+03]" ~ "3" )) Prey_count_df_FE$Prey_weight_class<-as.numeric(Prey_count_df_FE$Prey_weight_class) Prey_countLmmodel_FE <- lme(log10_Prey_count ~ Avg_SST + Mean_effort + Avg_SST*Mean_effort + Avg_SST*Prey_weight_class + Mean_effort*Prey_weight_class, random = list(ICES = ~ 1, year_since = ~1, Quarter = ~1, chl = ~1, meanSal = ~1, Depth = ~1, Pred_taxa = ~1 ), data = Prey_count_df_FE, na.action = na.exclude) Prey_countLmmodel_FE_poly <- lme(I(log10_Prey_count^2) ~ I(Avg_SST^2) + I(Mean_effort^2) + I(Avg_SST^2)*I(Mean_effort^2)+ I(Avg_SST^2)*I(Prey_weight_class^2) + I(Mean_effort^2)*I(Prey_weight_class^2), random = list(ICES = ~ 1, year_since = ~1, Quarter = ~1, chl = ~1, meanSal = ~1, Depth = ~1, Pred_taxa = ~1 ), data = Prey_count_df_FE, na.action = na.exclude) AIC(Prey_countLmmodel_FE, Prey_countLmmodel_FE_poly) #Model Results summary(Prey_countLmmodel_FE) r.squaredGLMM(Prey_countLmmodel_FE) setwd("C:/Users/as21769/OneDrive - University of Essex/Documents/Essex/Phd/2023/PPMR chapter/Resubmission April 2025") saveRDS(Prey_countLmmodel_FE, "Prey_countLmmodel_FE.rdata") Prey_countLmmodel_FE<-readRDS("Prey_countLmmodel_FE.rdata") #Model Predictions new_data_count <- expand.grid(Avg_SST = SST_quantiles, Mean_effort = FE_medians, Prey_weight_class = c(1, 3)) predictions_countInt<-predict(Prey_countLmmodel_FE, newdata = new_data_count, level = 0, se.fit = TRUE, type = "link") %>% as.data.frame() dataInt_count<-as.data.frame(cbind(new_data_count, predictions_countInt)) dataInt_count$Mean_effort<-as.factor(dataInt_count$Mean_effort) dataInt_count$Avg_SST<-as.numeric(dataInt_count$Avg_SST) dataInt_count$fit<-as.numeric(dataInt_count$fit) dataInt_count$Prey_weight_class<-as.factor(dataInt_count$Prey_weight_class) dataInt_count<- dataInt_count %>% mutate(FE_class = case_when(Mean_effort == FE_medians[1] ~ "1", Mean_effort == FE_medians[2] ~ "2", Mean_effort== FE_medians[3] ~ "3" )) p_Prey_count_Int <- ggplot() + geom_point(data = Prey_count_df_FE, aes(x = Avg_SST, y = log10_Prey_count, color = as.factor(FE_class)), alpha = 0.25, size = 1) + geom_smooth(method = "lm", aes(x = Avg_SST, y = fit, color = as.factor(FE_class), linetype = Prey_weight_class), data = subset(dataInt_count, FE_class %in% c("1", "3")), size = 1.5, se = TRUE) + scale_color_manual(labels = c('Low', 'Medium', 'High'), values = c("#76c68f", "grey", "#bd7ebe")) + labs( x = "Temperature (°C)", y = expression(bold(Log[10] ~ Prey ~ Count)), color = "Mean Fishing Effort" ) + theme_classic() + theme( axis.text.x = element_text(size = 24, colour = "black"), axis.text.y = element_text(size = 24, colour = "black"), axis.title = element_text(face = "bold", size = 25), plot.subtitle = element_text(size = 26, face = "bold", family = "arial") )+scale_linetype_manual(labels = c("Small", "Large"), values = c("solid", "dashed"))+ theme(legend.position = "none")+ labs(subtitle = "c", size = 10, fontface = "bold")#+ # stat_regline_equation( # data = subset(dataInt_count, FE_class %in% c("1", "3")), # aes(x = Avg_SST, y = fit, label = ..eq.label.., color =FE_class, linetype =Prey_weight_class), # label.y = c(3.6, 3.9, 4.2, 4.5), # label.x = 12.2 # ) p_Prey_count_Int setwd("C:/Users/as21769/OneDrive - University of Essex/Documents/Essex/Phd/2023/PPMR chapter/FINAL PPMR PUBLICATION") model_data <- eval(Prey_countLmmodel_FE$call$data, environment(formula(Prey_countLmmodel_FE))) mf <- model.frame(formula(Prey_countLmmodel_FE), data = model_data, na.action = na.omit) mf<-unique(na.omit(mf)) write.csv(mf, "Prey_countLmmodel_FE.csv") plot_data <- ggplot_build(p_Prey_count_Int)$data[[2]] write.csv(plot_data, "Figure 1 Prey_countLmmodel_FE.csv") ##Prey Richness df_cleaned_PPMR$Species<-1 Prey_Rich_All_Data<-df_cleaned_PPMR %>% group_by(stomach_id,Prey_taxa) %>% dplyr::summarize(Prey_count= sum(Prey_count), Prey_weight = mean(Prey_weight), Avg_SST = mean(Avg_SST), meanSal = mean(meanSal), chl = mean(chl), year_since = year_since, ICES = ICES, Depth = mean(Depth), Quarter = Quarter, Richness = sum(Species), Pred_taxa = Pred_taxa) Prey_Rich_All_Data$Prey_weight_class = as.numeric(cut_number(Prey_Rich_All_Data$Prey_weight, 3)) Prey_Rich_All_Data$log10_Richness<-log10(Prey_Rich_All_Data$Richness) RichnessLmmodel <- lme(log10_Richness ~ Avg_SST + Avg_SST*Prey_weight_class, random = list(ICES = ~ 1, year_since = ~1, Quarter = ~1, chl = ~1, meanSal = ~1, Depth = ~1, Pred_taxa = ~1 ), data = Prey_Rich_All_Data, na.action = na.exclude) RichnessLmmodel_poly <- lme(log10_Richness ~ I(Avg_SST^2) + I(Avg_SST^2)*I(Prey_weight_class^2), random = list(ICES = ~ 1, year_since = ~1, Quarter = ~1, chl = ~1, meanSal = ~1, Depth = ~1, Pred_taxa = ~1 ), data = Prey_Rich_All_Data, na.action = na.exclude) AIC(RichnessLmmodel, RichnessLmmodel_poly) setwd("C:/Users/as21769/OneDrive - University of Essex/Documents/Essex/Phd/2023/PPMR chapter/Resubmission April 2025") saveRDS(RichnessLmmodel_poly, "RichnessLmmodel_poly.rdata") RichnessLmmodel_poly<-readRDS("RichnessLmmodel_poly.rdata") #Model Results summary(RichnessLmmodel_poly) r.squaredGLMM(RichnessLmmodel) #Model Predictions predictions_Rich<-predict(RichnessLmmodel_poly, newdata = Prey_Rich_All_Data, level = 0, se.fit = TRUE, type = "link") %>% as.data.frame() data_Rich<-cbind(Prey_Rich_All_Data, predictions_Rich) data_Rich$Avg_SST<-as.numeric(data_Rich$Avg_SST) data_Rich$log10_Richness<-as.numeric(data_Rich$log10_Richness) data_Rich$Prey_weight_class<-as.factor(data_Rich$Prey_weight_class) y.expression.count <- expression(bold(Log[10] ~Prey ~ Richness)) Prey_Rich_All_Data$Prey_weight_class<-as.factor(Prey_Rich_All_Data$Prey_weight_class) Prey_Rich_plot<-ggplot() + geom_point(aes(y=Prey_Rich_All_Data$log10_Richness, x= Prey_Rich_All_Data$Avg_SST,colour=Prey_Rich_All_Data$Prey_weight_class),group=Prey_Rich_All_Data$Prey_weight_class, alpha = 0.35, size = 1) + geom_smooth(method = "lm", formula = y ~ I(x^2), aes(x = Avg_SST, y = fit, color = as.factor(Prey_weight_class)), data = data_Rich, size = 1.5, se = TRUE) + scale_color_manual(labels=c('Small', 'Medium', 'Large'), values = pal3) + scale_fill_manual(values = pal3) + labs(x= "Temperature (°C)", y=y.expression.count, color="Prey\nweight class", fill="Prey weight\nclass") + theme_classic() + theme(axis.text.x=element_text(size=24,colour="black")) + theme(axis.text.y=element_text(size=24,colour="black")) + theme(axis.title=element_text(face="bold", size=25)) + theme(legend.title =element_text(face="bold", size=18)) + theme(legend.text =element_text(size=18))+ theme(legend.position = "none")+ labs(subtitle = "d", size=10,fontface = "bold") + theme(plot.subtitle = element_text(size = 26, face = "bold",family="arial")) #+ # stat_regline_equation( # data = data_Rich, # aes(x = Avg_SST, y = fit, label = ..eq.label.., color =Prey_weight_class), # label.y = c(3.6, 3.9, 4.2, 4.5), # label.x = 12.2, # formula = y ~ poly(x, 2, raw = TRUE) # ) Prey_Rich_plot setwd("C:/Users/as21769/OneDrive - University of Essex/Documents/Essex/Phd/2023/PPMR chapter/FINAL PPMR PUBLICATION") model_data <- eval(RichnessLmmodel_poly$call$data, environment(formula(RichnessLmmodel_poly))) mf <- model.frame(formula(RichnessLmmodel_poly), data = model_data, na.action = na.omit) mf<-unique(na.omit(mf)) write.csv(mf, "RichnessLmmodel_poly.csv") plot_data <- ggplot_build(Prey_Rich_plot)$data[[2]] write.csv(plot_data, "Figure 1 RichnessLmmodel_poly.csv") #Richness Interaction df_cleaned_PPMR_FE$Species<-1 Prey_Rich_FE_Data<-df_cleaned_PPMR_FE %>% group_by(stomach_id,Prey_taxa) %>% dplyr::summarize(Prey_count= sum(Prey_count), Prey_weight = mean(Prey_weight), Avg_SST = mean(Avg_SST), meanSal = mean(meanSal), chl = mean(chl), year_since = year_since, ICES = ICES, Depth = mean(Depth), Quarter = Quarter, Richness = sum(Species), Mean_effort = mean(Mean_effort), Pred_taxa = Pred_taxa) Prey_Rich_FE_Data$Prey_weight_class = as.numeric(cut_number(Prey_Rich_FE_Data$Prey_weight, 3)) Prey_Rich_FE_Data$log10_Richness<-log10(Prey_Rich_FE_Data$Richness) Prey_Rich_FE_Data$FE_class<-as.numeric(cut_number(Prey_Rich_FE_Data$Mean_effort,3)) RichnessLmmodel_FE <- lme(log10_Richness ~ Avg_SST + Mean_effort + Avg_SST*Mean_effort + Avg_SST*Prey_weight_class + Mean_effort*Prey_weight_class, random = list(ICES = ~ 1, year_since = ~1, Quarter = ~1, chl = ~1, meanSal = ~1, Depth = ~1, Pred_taxa = ~1 ), data = Prey_Rich_FE_Data, na.action = na.exclude) RichnessLmmodel_FE_poly <- lme(log10_Richness ~ I(Avg_SST^2) + I(Mean_effort^2) + I(Avg_SST^2)*I(Mean_effort^2) + I(Avg_SST^2)*I(Prey_weight_class^2) + I(Mean_effort^2)*I(Prey_weight_class^2), random = list(ICES = ~ 1, year_since = ~1, Quarter = ~1, chl = ~1, meanSal = ~1, Depth = ~1, Pred_taxa = ~1 ), data = Prey_Rich_FE_Data, na.action = na.exclude) AIC(RichnessLmmodel_FE, RichnessLmmodel_FE_poly) setwd("C:/Users/as21769/OneDrive - University of Essex/Documents/Essex/Phd/2023/PPMR chapter/Resubmission April 2025") saveRDS(RichnessLmmodel_FE_poly, "RichnessLmmodel_FE_poly.rdata") RichnessLmmodel_FE_poly<-readRDS("RichnessLmmodel_FE_poly.rdata") #Model Results summary(RichnessLmmodel_FE_poly) r.squaredGLMM(RichnessLmmodel_FE_poly) coef(RichnessLmmodel_FE_poly) #Model Predictions new_data_Rich <- expand.grid(Avg_SST = SST_quantiles, Mean_effort = FE_medians, Prey_weight_class = c(1, 3)) predictions_RichInt<-predict(RichnessLmmodel_FE_poly, newdata = new_data_Rich, level = 0, se.fit = TRUE, type = "link") %>% as.data.frame() dataInt_Rich<-as.data.frame(cbind(new_data_Rich, predictions_RichInt)) dataInt_Rich$Mean_effort<-as.factor(dataInt_Rich$Mean_effort) dataInt_Rich$Avg_SST<-as.numeric(dataInt_Rich$Avg_SST) dataInt_Rich$fit<-as.numeric(dataInt_Rich$fit) dataInt_Rich$Prey_weight_class<-as.factor(dataInt_Rich$Prey_weight_class) dataInt_Rich<- dataInt_Rich %>% mutate(FE_class = case_when(Mean_effort == FE_medians[1] ~ "1", Mean_effort == FE_medians[2] ~ "2", Mean_effort== FE_medians[3] ~ "3" )) p_Prey_Rich_Int <- ggplot() + geom_point(data = Prey_Rich_FE_Data, aes(x = Avg_SST, y = log10_Richness, color = as.factor(FE_class)), alpha = 0.25, size = 1) + geom_smooth(method = "lm", formula = y ~ I(x^2), aes(x = Avg_SST, y = fit, color = as.factor(FE_class), linetype = Prey_weight_class), data = subset(dataInt_Rich, FE_class %in% c("1", "3")), size = 1.5, se = TRUE) + scale_color_manual(labels = c('Low', 'Medium', 'High'), values = c("#76c68f", "grey", "#bd7ebe")) + labs( x = "Temperature (°C)", y = expression(bold(Log[10] ~ Prey ~ Richness)), color = "Mean Fishing Effort" ) + theme_classic() + theme( axis.text.x = element_text(size = 24, colour = "black"), axis.text.y = element_text(size = 24, colour = "black"), axis.title = element_text(face = "bold", size = 25), plot.subtitle = element_text(size = 26, face = "bold", family = "arial") )+scale_linetype_manual(labels = c("Small", "Large"), values = c("solid", "dashed"))+ labs(subtitle = "d", size = 10, fontface = "bold")#+ theme(legend.position = "none")#+ # stat_regline_equation( # data = subset(dataInt_Rich, FE_class %in% c("1", "3")), # aes(x = Avg_SST, y = fit, label = ..eq.label.., color =FE_class, linetype =Prey_weight_class), # label.y = c(3.6, 3.9, 4.2, 4.5), # label.x = 12.2, # formula = y ~ poly(x, 2, raw = TRUE) #) p_Prey_Rich_Int setwd("C:/Users/as21769/OneDrive - University of Essex/Documents/Essex/Phd/2023/PPMR chapter/FINAL PPMR PUBLICATION") model_data <- eval(RichnessLmmodel_FE_poly$call$data, environment(formula(RichnessLmmodel_FE_poly))) mf <- model.frame(formula(RichnessLmmodel_FE_poly), data = model_data, na.action = na.omit) mf<-unique(na.omit(mf)) write.csv(mf, "RichnessLmmodel_FE_poly.csv") plot_data <- ggplot_build(p_Prey_Rich_Int)$data[[2]] write.csv(plot_data, "Figure 1 RichnessLmmodel_FE_poly.csv") ## All nmds # Clean and prepare NMDS input data NMDSDf <- df_cleaned_PPMR %>% transmute( Stomach = as.character(haul_id), Prey_taxa, Prey_count = as.numeric(Prey_count), Prey_weight = as.numeric(Prey_weight), Avg_SST = as.numeric(Avg_SST) ) %>% mutate(Prey_count = ifelse(is.na(Prey_count), 1, Prey_count)) %>% na.omit() NMDSDf$Stomach <- as.factor(NMDSDf$Stomach) data <- NMDSDf %>% group_by(Stomach, Prey_taxa) %>% dplyr::summarize( Prey_count = sum(Prey_count, na.rm = TRUE), Prey_weight = mean(Prey_weight, na.rm = TRUE), Avg_SST = mean(Avg_SST, na.rm = TRUE), .groups = "drop" ) %>% mutate(Prey_taxa = str_replace(Prey_taxa, " ", "_")) ENVDf <- df_cleaned_PPMR %>% transmute( Stomach = as.character(haul_id), Avg_SST = as.numeric(Avg_SST), Prey_Body_Mass = as.numeric(Prey_weight) ) %>% group_by(Stomach) %>% dplyr::summarize( Avg_SST = mean(Avg_SST, na.rm = TRUE), Prey_Body_Mass = mean(Prey_Body_Mass, na.rm = TRUE), .groups = "drop" ) stomachs <- data.frame(Stomach = unique(data$Stomach)) env_df <- merge(stomachs, ENVDf, by = "Stomach", all.x = TRUE) rownames(env_df) <- env_df$Stomach # Compute prey weight classes sum_data <- data %>% group_by(Prey_taxa) %>% dplyr::summarize(Prey_weight = mean(Prey_weight, na.rm = TRUE), .groups = "drop") %>% mutate(Prey_weight_class = as.numeric(cut_number(Prey_weight, 3))) # Reformat to wide matrix: Stomach x Prey_taxa data_matrix <- data %>% group_by(Stomach, Prey_taxa) %>% dplyr::summarize(Prey_count = sum(Prey_count, na.rm = TRUE), .groups = "drop") %>% acast(Stomach ~ Prey_taxa, value.var = "Prey_count", fun.aggregate = mean) data_matrix[is.na(data_matrix)] <- 0 data_matrix <- data_matrix[rowSums(data_matrix) > 0, colSums(data_matrix) > 0] # Filter env_df to match matrix rows env_df <- env_df[rownames(data_matrix), ] # Recheck alignment stopifnot(nrow(data_matrix) == nrow(env_df)) trymax = 999 # up to 999 if not converging dims = 2 # number of dimensions (I like to keep this down else results become long and hard to interpret) set.seed(5) # means results reproducible nmds = metaMDS(data_matrix, dist='euclidian', binary=F, k = dims, sratmax =0.9999999, trymax=12, maxit=999) setwd("C:/Users/as21769/OneDrive - University of Essex/Documents/Essex/Phd/2023/PPMR chapter/Resubmission April 2025") saveRDS(nmds, "nmds.rdata") nmds<-readRDS("nmds.rdata") nfit(nmds) vegdist(data1) species_score = as.data.frame(scores(nmds)$species) data.scores_sites = as.data.frame(scores(nmds)$sites) en_coord_cont = as.data.frame(scores(envfit, "vectors")) sum_data<-subset(sum_data, Prey_taxa %in% rownames(species_score)) envfit<-envfit(nmds ~ Prey_Body_Mass +Avg_SST, data = env_df) spp.scrs <- as.data.frame(scores(envfit, display = "vectors")) sum_data$Prey_weight_class<-as.factor(sum_data$Prey_weight_class) SST_mds<-ggplot() + geom_point(data=species_score,aes(x=NMDS1,y=NMDS2,colour=sum_data$Prey_weight_class,size=sum_data$Prey_weight_class))+ coord_equal() + labs(color="Prey\nweight class")+ theme_classic()+ geom_segment(data = spp.scrs, aes(x = 0, xend = NMDS1, y = 0, yend = NMDS2), arrow = arrow(length = unit(0.25, "cm")), colour = "black", size = 1)+ geom_text(data = spp.scrs, aes(x = NMDS1, y = NMDS2, label = c("Prey Body Mass", "SST")), size = 2.5,fontface = "bold", hjust=c("centre"), vjust = c(-0.5, 2.8))+ theme(axis.text=element_text(size=14,colour="black")) + theme(axis.title=element_text(face="bold", size=25)) + theme(legend.title =element_text(face="bold", size=24)) + theme(legend.text =element_text(size=22)) + theme(plot.subtitle = element_text(size = 26, face = "bold",family="arial")) + scale_colour_manual(labels=c('Small', 'Medium', 'Large'), values=pal3, drop = F, name = "Prey\nweight class") + scale_size_discrete(labels=c('Small', 'Medium', 'Large'), drop = F, name = "Prey\nweight class")+ labs(subtitle = "e", size=10,fontface = "bold")+ scale_colour_manual(labels=c('Small', 'Medium', 'Large'), values=pal3, drop = F, name = "Prey\nweight class")+ theme(legend.position = "none") SST_mds ##Fishing Effort NMDS NMDSDf_FE <- df_cleaned_PPMR_FE %>% transmute( Stomach = as.character(haul_id), Prey_taxa, Prey_count = as.numeric(Prey_count), Prey_weight = as.numeric(Prey_weight), Avg_SST = as.numeric(Avg_SST) ) %>% replace_na(list(Prey_count = 1)) %>% na.omit() NMDSDf_FE$Stomach <- as.factor(NMDSDf_FE$Stomach) data_FE <- NMDSDf_FE %>% group_by(Stomach, Prey_taxa) %>% dplyr::summarize( Prey_count = sum(Prey_count, na.rm = TRUE), Prey_weight = mean(Prey_weight, na.rm = TRUE), Avg_SST = mean(Avg_SST, na.rm = TRUE), .groups = "drop" ) %>% mutate(Prey_taxa = str_replace_all(Prey_taxa, " ", "_")) # Environmental data (SST, body mass, fishing effort) ENVDf_FE <- df_cleaned_PPMR_FE %>% transmute( Stomach = as.character(haul_id), Avg_SST = as.numeric(Avg_SST), Prey_Body_Mass = as.numeric(Prey_weight), Mean_effort = as.numeric(Mean_effort) ) %>% group_by(Stomach) %>% dplyr::summarize( Avg_SST = mean(Avg_SST, na.rm = TRUE), Prey_Body_Mass = mean(Prey_Body_Mass, na.rm = TRUE), Mean_effort = mean(Mean_effort, na.rm = TRUE), .groups = "drop" ) # Merge env data with feeding individuals env_df_FE <- data_FE %>% distinct(Stomach) %>% left_join(ENVDf_FE, by = "Stomach") %>% column_to_rownames("Stomach") # Prey weight classes sum_data_FE <- data_FE %>% group_by(Prey_taxa) %>% dplyr::summarize(Prey_weight = mean(Prey_weight, na.rm = TRUE), .groups = "drop") %>% mutate(Prey_weight_class = as.factor(cut_number(Prey_weight, 3))) # Prepare wide data for NMDS data1_FE <- data_FE %>% group_by(Stomach, Prey_taxa) %>% dplyr::summarize(Prey_count = sum(Prey_count, na.rm = TRUE), .groups = "drop") %>% pivot_wider(names_from = Prey_taxa, values_from = Prey_count, values_fill = 0) %>% column_to_rownames("Stomach") # Remove empty rows/cols data1_FE <- data1_FE[rowSums(data1_FE) > 0, colSums(data1_FE) > 0] env_df_FE <- env_df_FE[rownames(env_df_FE) %in% rownames(data1_FE), ] # Run NMDS set.seed(5) nmds_FE <- metaMDS(data1_FE, dist = "euclidean", k = 2, trymax = 12, maxit = 999) setwd("C:/Users/as21769/OneDrive - University of Essex/Documents/Essex/Phd/2023/PPMR chapter/Resubmission April 2025") saveRDS(nmds_FE, "nmds_FE.rdata") nmds_FE<-readRDS("nmds_FE.rdata") # Scores species_score_FE <- as.data.frame(scores(nmds_FE, "species")) data_scores_sites_FE <- as.data.frame(scores(nmds_FE, "sites")) # Filter prey taxa for species scores sum_data_FE <- sum_data_FE %>% filter(Prey_taxa %in% rownames(species_score_FE)) # Environmental fitting envfit_FE <- envfit(nmds_FE ~ Avg_SST + Prey_Body_Mass + Mean_effort, data = env_df_FE) spp_scrs_FE <- as.data.frame(scores(envfit_FE, display = "vectors")) # Plot SST_mds_FE <- ggplot() + geom_point(data = species_score_FE, aes(x = NMDS1, y = NMDS2, colour = sum_data_FE$Prey_weight_class, size = sum_data_FE$Prey_weight_class)) + coord_equal() + scale_color_manual(labels = c("Small", "Medium", "Large"), values = pal3, name = "Prey\nweight class", drop = FALSE) + scale_size_discrete(labels = c("Small", "Medium", "Large"), name = "Prey\nweight class", drop = FALSE) + labs(color = "Prey\nweight class") + geom_segment(data = spp_scrs_FE, aes(x = 0, xend = NMDS1, y = 0, yend = NMDS2), arrow = arrow(length = unit(0.25, "cm")), colour = "black", size = 1) + geom_text(data = spp_scrs_FE, aes(x = NMDS1, y = NMDS2, label = c(" Temperature", " Prey Body Mass", " Fishing effort")), size = 2.5, fontface = "bold", hjust = c("bottom", "bottom", "top")) + theme_classic() + theme( axis.text = element_text(size = 24, colour = "black"), axis.title = element_text(face = "bold", size = 25), legend.title = element_text(face = "bold", size = 24), legend.text = element_text(size = 22), plot.subtitle = element_text(size = 26, face = "bold", family = "arial") ) + labs(subtitle = "e")#+ theme(legend.position = "none") SST_mds_FE