--- title: "Deepwater Red Sea" author: "Andrew Temple" date: "2026-02-08" output: html_document editor_options: chunk_output_type: console --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE) ``` # Set Up Load libraries ```{r} library(tidyverse) library(cowplot) library(ggridges) library(rfishbase) library(worrms) library(robis) library(sf) library(rnaturalearth) library(geosphere) library(mice) library(howManyImputations) library(VIM) library(naniar) library(ape) library(coRanking) library(ggtree) library(brms) library(tidybayes) library(loo) library(corrplot) library(wesanderson) library(magick) library(cluster) library(phytools) library(GGally) library(grid) library(vegan) library(geometry) ``` # Data Compilation ## Species compilation Get species list from Fishbase for ALL deep-water species globally. we define deep water as those >200m depth (e.g. beyond the shelf). ```{r} sp_all <- load_taxa() %>% left_join( estimate() %>% dplyr::select(c("SpecCode", "DepthMin", "DepthMax", "MaxLengthTL", "Troph", "K", "FeedingPath")) ) %>% left_join( species() %>% dplyr::select(c("SpecCode", "BodyShapeI", "DemersPelag", "Saltwater", "Brack")) ) %>% left_join( reproduction() %>% dplyr::select(c("SpecCode", "Fertilization")) ) %>% ## depth filter filter(DepthMax > 200) %>% ## brackish and/or saltwater filter filter(Brack == "1" | Saltwater == "1") %>% ## remove columns we no longer need dplyr::select(-c(Brack, Saltwater)) ## 6354 species ``` Get their WORMS AlphaIDs, will be used to extract OBIS data ```{r} get_worms_id <- function(name) { tryCatch({ id <- wm_name2id(name = name) if (is.null(id)) { return(NA) } else { return(id) } }, error = function(e) { # Handle cases where no match is found or an error occurs return(NA) }) } sp_all$AphiaID <- sapply(sp_all$Species, get_worms_id) ## Filter to only species with AphiaID (WORMS) sp_all <- sp_all %>% drop_na(AphiaID) ``` Load fishbase list for RS/AS_SC and append it to the deep sea list ```{r} sp <- read.csv("Species_list.csv") sp_all <- sp_all %>% left_join(sp) %>% mutate(across(all_of(c("Red_Sea", "AS_SC")), ~replace_na(.x, 0))) ## clean rm(sp) ``` ## OBIS Distribution Data Get occurrences for all the species listed on OBIS ```{r} ## OBIS: EXTRACTION DONE 31 MAY 2026 records_obis <- lapply(sp_all$AphiaID, function(sp) { message("Downloading: ", sp) tryCatch( robis::occurrence(taxonid = sp) %>% mutate(AphiaID = sp) %>% select(AphiaID, decimalLongitude, decimalLatitude), error = function(e) NULL ) }) records_obis <- bind_rows(records_obis) ``` Clean OBIS data and also subset it for indian ocean species ```{r} ## remove Red Sea records from the data set records_obis <- records_obis %>% filter( !(decimalLongitude >= 32 & decimalLongitude <= 43.5 & decimalLatitude >= 12.54 & decimalLatitude <= 30) ) ## clean OBIS records to remove land points land <- ne_countries(scale = "medium", returnclass = "sf") records_obis <- st_as_sf(records_obis, coords = c("decimalLongitude", "decimalLatitude"), crs = 4326, remove = FALSE) on_land <- st_intersects(records_obis, land) records_obis <- records_obis[lengths(on_land) == 0, ] ## get an Indian Ocean only subset indian_ocean <- st_read("GOaS_v1_20211214/goas_v01.shp") %>% filter(name == "Indian Ocean") io <- st_intersects(records_obis, indian_ocean) records_io <- records_obis[lengths(io) == 1, ] ## convert back to dataframes records_obis <- records_obis %>% st_drop_geometry() records_io <- records_io %>% st_drop_geometry() ## clean up rm(on_land, io, indian_ocean, land) gc() ``` Calculate min and max latitude plus the total number of sightings/reports (as a proxy for discoverability) ```{r} ## min/max lat lat_ranges <- records_obis %>% group_by(AphiaID) %>% summarise( min_lat = min(decimalLatitude, na.rm = TRUE), max_lat = max(decimalLatitude, na.rm = TRUE), sightings = n(), .groups = "drop" ) sp_all <- sp_all %>% left_join(lat_ranges) ## if sightings = NA, set to 0 - its not missign data, it is true absence (on OBIS) sp_all <- sp_all %>% mutate(sightings = ifelse(is.na(sightings), 0, sightings)) ## clean rm(lat_ranges) ``` Calculate the nearest point to the Bab-el-Mendab straight (based on Indian Ocean points alone) ```{r} ## extract min and max latitude outside of Red Sea. we set the Bab al-Mandab Strait as = 12.54, 43.50 bab_lat <- 12.54 bab_lon <- 43.50 bab_dist <- records_io %>% mutate( dist_bab_m = distHaversine( cbind(decimalLongitude, decimalLatitude), c(bab_lon, bab_lat) ), dist_bab_km = dist_bab_m / 1000 ) bab_dist <- bab_dist %>% group_by(AphiaID) %>% summarise( min_dist_bab = min(dist_bab_km, na.rm = TRUE), .groups = "drop" ) ## join to data frame sp_all <- sp_all %>% left_join(bab_dist) ## For Red Se endemics set distance to be 1. sp_all <- sp_all %>% mutate(min_dist_bab = ifelse(sightings == 0 & is.na(min_dist_bab) & Red_Sea == 1 & Species != "Arnoglossus arabicus" & Species != "Saurenchelys lateromaculata", 1, min_dist_bab)) ## clean rm(bab_lat, bab_lon, bab_dist, records_obis, records_io) gc() ``` ## clean up data set Before we impute data for the analysis, we need to finalise the data Ensure basin assignments are all correct ```{r} sp_all <- sp_all %>% ## IO must include all Arabian Seas, Somali Current, and Red Sea species, along with any having OBIS points there mutate(IO = ifelse(!is.na(min_dist_bab) | AS_SC == 1 | Red_Sea == 1, 1, 0)) ``` we first clean up the data and calculate derived variables. ```{r} # make dataset that will become final, update classes that need to be updated first sapply(sp_all, class) sp_all <- sp_all %>% mutate(Fertilization = as.factor(ifelse(Fertilization == "internal (oviduct)", "internal", ifelse(Fertilization == "external", "external", NA)))) %>% mutate(Phylum = "Chordata") %>% dplyr::select(c(SpecCode, AphiaID, Phylum, SuperClass, Class, Order, Family, Genus, Species, ## Phylogeny IO, AS_SC, Red_Sea, sightings, ## Presence min_lat, max_lat, DepthMin, DepthMax, ## Distribution min_dist_bab, ## Access to RS MaxLengthTL, Troph, K, BodyShapeI, DemersPelag, Fertilization, FeedingPath ## Ecology )) ## clean up categories which are messy unique(sp_all$BodyShapeI) unique(sp_all$DemersPelag) sp_all <- sp_all %>% mutate(BodyShapeI = ifelse(BodyShapeI == "Elongated", "elongated", BodyShapeI), DemersPelag = case_when( DemersPelag %in% c("bathydemersal", "demersal", "reef-associated") ~ "demersal", DemersPelag %in% c("benthopelagic") ~ "benthopelagic", DemersPelag %in% c("bathypelagic", "pelagic-oceanic", "pelagic-neritic", "pelagic") ~ "pelagic", TRUE ~ NA_character_ )) ``` Calculate derived variables that need to be imputed. We will impute median latitude and latitude range because min/max latitude are related and jointly missing. Imputing min/max will be a mess and the derived variables from imputations will be deeply unreliable. median and range are more dissociated from one another so will be less problematic. depth range and median depth can also be derived now, and we will derive them from imputations again later (because max depth isn't missing which will remove the issues in min depth imputation.) ```{r} #t <- sp_all %>% mutate(x = DepthMax - DepthMin) ## one species has the min and max inverted... ## fix this and re-calculate! sp_all <- sp_all %>% mutate(DepthMin = ifelse(Species == "Pogonophryne sarmentifera", 1036, DepthMin), DepthMax = ifelse(Species == "Pogonophryne sarmentifera", 1157, DepthMax),) ## recalculate sp_all <- sp_all %>% mutate( med_depth = ((DepthMax - DepthMin)/2) + DepthMin, depth_range = DepthMax - DepthMin, med_lat = (min_lat + max_lat)/2, lat_range = max_lat - min_lat ) ``` Save the data ```{r} #write.csv(sp_all, "final_species_data.csv", row.names = FALSE) ``` # Analysis: Imputation and Phylogeny Filter records as needed ```{r} sp_all <- read.csv("final_species_data.csv") ## ensure classes are correct sapply(sp_all, class) sp_all <- sp_all %>% mutate(across(where(is.character), as.factor)) %>% mutate(across(where(is.integer), as.numeric)) ``` We need to transform data before imputation, since imputation uses linear models highly skewed variables are not ideal ```{r} ## check variables to be included in the discovery and functional space models only hist(sp_all$sightings) ## skew hist(sp_all$min_lat) hist(sp_all$max_lat) hist(sp_all$med_lat) hist(sp_all$DepthMin) ## skew hist(sp_all$DepthMax) ## skew hist(sp_all$min_dist_bab) hist(sp_all$MaxLengthTL) ## skew hist(sp_all$Troph) hist(sp_all$K) ## skew hist(sp_all$med_depth) ## skew hist(sp_all$depth_range) ## skew hist(sp_all$lat_range) ## skew ## log transform those where it is needed, we can backtransform later ## log10(var + 1) all just for consistency (some have 0s) sp_all <- sp_all %>% mutate(sightings = log10(sightings+1), DepthMin = log10(DepthMin+1), DepthMax = log10(DepthMax+1), MaxLengthTL = log10(MaxLengthTL+1), K = log10(K+1), med_depth = log10(med_depth+1), depth_range = log10(depth_range+1), lat_range = log10(lat_range+1) ) ``` filter to only columns we want to take forwards ```{r} sp_all <- sp_all %>% dplyr::select(IO, AS_SC, Red_Sea, Phylum, Class, Order, Family, Genus, Species, sightings, DepthMin, DepthMax, min_dist_bab, MaxLengthTL, Troph, K, BodyShapeI, DemersPelag, Fertilization, FeedingPath, med_depth, depth_range, med_lat, lat_range) ``` ## imputation Check missingness patterns ```{r} # Pattern and proportion overview only for variables we will impute/use in imputation md.pattern(sp_all %>% dplyr::select(-c(Phylum, Order, Family, Genus, Species, depth_range, med_depth)), rotate.names = TRUE) aggr(sp_all %>% dplyr::select(-c(Phylum, Order, Family, Genus, Species, depth_range, med_depth)), col = c("steelblue","red"), numbers = TRUE, sortVars = TRUE) # Little's MCAR test (H0: MCAR) mcar_test(sp_all %>% dplyr::select(-c(Phylum, Order, Family, Genus, Species, depth_range, med_depth))) # p < 0.05 = reject MCAR ``` Set methods for imputation ```{r} ini <- mice(sp_all, maxit = 0) ini$method ## Fully control what gets imputed meth <- ini$method ## ensure all variables for imputation are correctly specified as pmm/logreg meth ## med_depth and depth_range will be imputed passively. ## these need to account for the log transformations we already did! ## we have to revert the log transformation, get the value, then transform it back into log scale (because everything else is in log10(+1) scale) meth["depth_range"] <- "~I(log10((10^DepthMax - 10^DepthMin) + 1))" meth["med_depth"] <- "~I(log10(((10^DepthMin + 10^DepthMax - 2)/2) + 1))" ## DepthMin needs a custom imputation becasue DepthMin cannot ever be > DepthMax, standard pmm violates this. mice.impute.pmm.bounded <- function(y, ry, x, wy = NULL, ...) { # Standard PMM draw vals <- mice.impute.pmm(y, ry, x, wy, ...) # DepthMax for recipient rows (those needing imputation) depthmax_recipients <- x[!ry, "DepthMax"] # Replace any violation by resampling from valid observed donors violations <- vals > depthmax_recipients if (any(violations)) { for (i in which(violations)) { valid_donors <- y[ry & y < depthmax_recipients[i]] if (length(valid_donors) == 0) { # Fallback: use DepthMax - 1 if no valid donors exist vals[i] <- depthmax_recipients[i] - 1 } else { vals[i] <- sample(valid_donors, 1) } } } return(vals) } # Assign custom method meth["DepthMin"] <- "pmm.bounded" ## check meth ``` Set prediction matrix for imputation ```{r} ## Must include all key variables to be used in the later model (congeniality principle) pred <- ini$predictorMatrix pred ## don't allow derived variables to be predictors for their parents pred["DepthMin", "depth_range"] <- 0 pred["DepthMin", "med_depth"] <- 0 ## only allow Class as a taxonomic predictor pred[, c("Phylum", "Order", "Family", "Genus", "Species")] <- 0 pred ``` Test run imputations to estimate FMI (Fraction of Missing Information) and therefore calculate final number of imputations needed Use the von Hippel quadratic rule. https://statisticalhorizons.com/how-many-imputations/ ```{r} imp_pilot <- mice( sp_all, m = 20, ## base recommended by von Hippel 20 method = meth, predictorMatrix = pred, maxit = 25, ## it usually 25 seed = 2009 ) ## check imputations look ok imp_1 <- complete(imp_pilot, action = "all") imp_1 <- imp_1[[1]] # create a 'mask' of rows to exclude from analysis that will not be used in the discovery model #ignore <- imp_pilot$data$IO == 0 regional <- (imp_pilot$data$AS_SC == 1) ## get FMI ## keep in mind that DepthMin, med_depth, depth_range, lat_range, MaxLengthTL, sightings, and K are already log10 transformed fit_pilot <- pool( with( imp_pilot, glm(Red_Sea ~ scale(DepthMin) + scale(med_lat) + scale(min_dist_bab) + ## access scale(depth_range) + scale(lat_range) + ## tolerance scale(MaxLengthTL) + scale(Troph) + ## ecology BodyShapeI + DemersPelag + Fertilization + ## ecology scale(sightings), family = binomial, subset = regional # <-- only use rows for the species pool ) ) ) how_many_imputations(fit_pilot, cv = .05, alpha = .05) how_many_imputations(fit_pilot, cv = .01, alpha = .05) ## 50 imputations needed for .05 and 1213 needed for .01! ## clean up rm(imp_pilot, fit_pilot, regional) gc() ``` Now impute full set ```{r} d_imp <- mice( sp_all, m = 50, method = meth, predictorMatrix = pred, maxit = 25, seed = 2009 ) plot(d_imp) densityplot(d_imp, ~DepthMin) ## clean up rm(meth, pred, ini) ``` ## Phylogeny Create full phylogeny for the species pool based on taxonomy ```{r} sp_all <- sp_all %>% mutate(across(where(is.character), as.factor)) # Build tree from taxonomy tax_tree <- as.phylo.formula(~ Phylum/Class/Order/Family/Genus/Species, data = sp_all %>% filter(AS_SC == 1)) ## make tree ultametric (i.e. from class to species is equal length regardless of number of splits) tax_tree <- compute.brlen(tax_tree, method = "Grafen") plot(tax_tree) # Create phylogenetic covariance A <- ape::vcv(tax_tree) ``` # Analysis: Prediction Model ## Model Prepare data for analysis by creating full list (this is what brms accepts) and dropping non-Arabian Seas and non-Somali Current species ```{r} sp_final <- complete(d_imp, action = "all") ## We excluded species failing a priori biogeographic filters from the final analysis, but retained them during multiple imputation to preserve trait covariance structure. ## now retain only those in geographic regions relevant to the analysis (Arabian Seas and Somali Current) ## make all characters variables factors sp_final <- lapply(sp_final, function(df) { df %>% filter(AS_SC == 1) %>% mutate(across(where(is.character), as.factor)) }) ``` Back-transform the variables we log transformed for imputation. Alter reference levels for factors where needed. ```{r} sp_final <- lapply(sp_final, function(df) { df %>% mutate(sightings = 10^sightings -1, DepthMin = 10^DepthMin -1, DepthMax = 10^DepthMax -1, MaxLengthTL = 10^MaxLengthTL -1, K = 10^K -1, med_depth = 10^med_depth -1, depth_range = 10^depth_range -1, lat_range = 10^lat_range -1 ) %>% mutate(BodyShapeI = relevel(factor(BodyShapeI), ref = "eel-like"), DemersPelag = relevel(factor(DemersPelag), ref = "demersal") ) }) ``` Set weakly informative priors for the model ```{r} ## observed baseline is 170 species from 432, but we believe the reality is higher. ## Based on baseline default prior would be ~ logit(170/432), aka log((170/432)/(1-(170/432))) = -0.4325461 ## We hypothesise missing lineages so need prior that expects closer ratio, we need it to be weak because we don't know how much by. ## Right now the ratio is .39, so lets set a weak prior closer to .5. logit .5 = 0 #prior(normal(0, 1), class = "Intercept") ## classical weak prior for slopes (which could go in either direction), but enough to moderate extremes #prior(normal(0, 1), class = "b") ## prior for phylogeny, regularize so that the model weights traits over phylogeny unless the phylogeny effect is strongly supported. #prior(exponential(2), class = "sd", group = "species") ## priors priors <- c( prior(normal(0, 0.7), class = "b"), # traits prior(normal(0, 1), class = "Intercept"), # baseline (adjust as discussed) prior(exponential(2), class = "sd", group = "Species") # phylogeny ) ``` Model use brm for pooling multiple data sets. Scale all variables and log10 transform any that are highly skewed ```{r} ## test run on one data set # Extract the first imputed dataset (change index for any other iteration) imp_1 <- sp_final[[1]] # Remove iteration prefix from column names (e.g. "X1.DepthMin" → "DepthMin") names(imp_1) <- gsub("^X\\d+\\.", "", names(imp_1)) ## test m00 <- brm( Red_Sea ~ scale(log10(DepthMin+1)) + scale(med_lat) + scale(min_dist_bab) + ## access scale(log10(depth_range+1)) + scale(log10(lat_range+1)) + ## tolerance scale(log10(MaxLengthTL+1)) + scale(Troph) + ## ecology BodyShapeI + DemersPelag + Fertilization + ## ecology scale(log10(sightings+1)) + ## discoverability (1|gr(Species, cov = A)), ## phylogenetic effect chains = 4, cores = 4, iter = 6000, warmup = 4000, threads = threading(3), prior = priors, data = imp_1, data2 = list(A = A), family = bernoulli(link = "logit"), control = list(adapt_delta = 0.99, max_treedepth = 15), seed = 2009 ) summary(m00) ## model m01 <- brm_multiple( Red_Sea ~ scale(log10(DepthMin+1)) + scale(min_dist_bab) + scale(med_lat) + ## access scale(log10(depth_range+1)) + scale(log10(lat_range+1)) + ## tolerance scale(log10(MaxLengthTL+1)) + scale(Troph) + ## ecology BodyShapeI + DemersPelag + Fertilization + ## ecology scale(log10(sightings+1)) + ## discoverability (1|gr(Species, cov = A)), ## phylogenetic effect chains = 4, cores = 4, iter = 6000, warmup = 4000, threads = threading(3), prior = priors, data = sp_final, data2 = rep(list(list(A = A)), length(sp_final)), family = bernoulli(link = "logit"), control = list(adapt_delta = 0.99, max_treedepth = 15), seed = 2009 ) summary(m01) bayes_R2(m01) ## 0.2010362 ## need to remove some elements of the model so as not to blow up the storage required! ## the below removes the individual imputation models, which is fine because we never need them, we only use the pooled model. m01$fits <- NULL rm(imp_1, m00) gc() ``` Check co-linearity ```{r} post <- as_draws_df(m01) pars <- post[, grepl("^b_", names(post))] cor_post <- cor(pars) round(cor_post, 2) corrplot(cor_post, method = "color", tl.cex = 0.7, addCoef.col = 'grey') ## colinearity only really seen between levels of the same factor and lat_range & sightings. Not at extreme levels. rm(pars, post, cor_post) ``` Check model ```{r} ## Posterior predictive check - can the model reproduce the observed discovery process? pp_check(m01) pp_check(m01, type = "bars") # good ## Influence/leverage diagnostics - loo(m01) ``` Check likelihood of effects ```{r} # Extract the combined posterior draws as a data frame draws <- as_draws_df(m01) colnames(draws) # Calculate probability for a specific parameter (e.g., b_Intercept) print(mean(draws$b_scalelog10DepthMinP1 < 0)) print(mean(draws$b_scalemin_dist_bab < 0)) print(mean(draws$b_scalemed_lat > 0)) print(mean(draws$b_scalelog10depth_rangeP1 < 0)) print(mean(draws$b_scalelog10lat_rangeP1 > 0)) print(mean(draws$b_scalelog10MaxLengthTLP1 < 0)) print(mean(draws$b_scaleTroph > 0)) print(mean(draws$b_BodyShapeIelongated < 0)) print(mean(draws$b_BodyShapeIfusiformDnormal < 0)) print(mean(draws$b_BodyShapeIother > 0)) print(mean(draws$b_BodyShapeIshortandDordeep < 0)) print(mean(draws$b_DemersPelagbenthopelagic < 0)) print(mean(draws$b_DemersPelagpelagic < 0)) print(mean(draws$b_Fertilizationinternal < 0)) print(mean(draws$b_scalelog10sightingsP1 < 0)) rm(draws) ``` ## Plot Plot the effects ```{r} ## Distance to Bab-el-Mendab p01 <- conditional_effects(m01, effects = "min_dist_bab", re_formula = NA, spaghetti = FALSE)$min_dist_bab %>% ggplot(aes(x = min_dist_bab, y = estimate__)) + geom_ribbon(aes(ymin = lower__, ymax = upper__), colour = NA, fill = "#3A9AB2", alpha = .25) + geom_line(colour = "#3A9AB2", linewidth = 2) + scale_y_continuous(limits = c(0,1), breaks = c(0,1)) + labs(x = "Bab-el-Mendab (km)", y = "Probability") + theme_classic() + theme(text = element_text(size = 6), axis.text.x = element_text(angle = 30, hjust = .5, vjust = .75)) ## Minimum Depth p02 <- conditional_effects(m01, effects = "DepthMin", re_formula = NA, spaghetti = FALSE)$DepthMin %>% ggplot(aes(x = DepthMin, y = estimate__)) + geom_ribbon(aes(ymin = lower__, ymax = upper__), colour = NA, fill = "#72B2BF", alpha = .25) + geom_line(colour = "#72B2BF", linewidth = 2) + scale_y_continuous(limits = c(0,1), breaks = c(0,1)) + labs(x = "Minimum Depth (m)", y = "Probability") + theme_classic() + theme(text = element_text(size = 6), axis.text.x = element_text(angle = 30, hjust = .5, vjust = .75)) ## Depth Range p03 <- conditional_effects(m01, effects = "depth_range", re_formula = NA, spaghetti = FALSE)$depth_range %>% ggplot(aes(x = depth_range, y = estimate__)) + geom_ribbon(aes(ymin = lower__, ymax = upper__), colour = NA, fill = "#95BBB1", alpha = .25) + geom_line(colour = "#95BBB1", linewidth = 2) + scale_y_continuous(limits = c(0,1), breaks = c(0,1)) + labs(x = "Depth Range (m)", y = "Probability") + theme_classic() + theme(text = element_text(size = 6), axis.text.x = element_text(angle = 30, hjust = .5, vjust = .75)) ## Median Latitude p04 <- conditional_effects(m01, effects = "med_lat", re_formula = NA, spaghetti = FALSE)$med_lat %>% ggplot(aes(x = med_lat, y = estimate__)) + geom_ribbon(aes(ymin = lower__, ymax = upper__), colour = NA, fill = "#ADC397", alpha = .25) + geom_line(colour = "#ADC397", linewidth = 2) + scale_y_continuous(limits = c(0,1), breaks = c(0,1)) + labs(x = "Median Latitude (degrees)", y = "Probability") + theme_classic() + theme(text = element_text(size = 6), axis.text.x = element_text(angle = 30, hjust = .5, vjust = .75)) ## Latitude Range p05 <- conditional_effects(m01, effects = "lat_range", re_formula = NA, spaghetti = FALSE)$lat_range %>% ggplot(aes(x = lat_range, y = estimate__)) + geom_ribbon(aes(ymin = lower__, ymax = upper__), colour = NA, fill = "#CAC96A", alpha = .25) + geom_line(colour = "#CAC96A", linewidth = 2) + scale_y_continuous(limits = c(0,1), breaks = c(0,1)) + labs(x = "Latitude Range (degrees)", y = "Probability") + theme_classic() + theme(text = element_text(size = 6), axis.text.x = element_text(angle = 30, hjust = .5, vjust = .75)) ## Maximum Length p06 <- conditional_effects(m01, effects = "MaxLengthTL", re_formula = NA, spaghetti = FALSE)$MaxLengthTL %>% ggplot(aes(x = MaxLengthTL, y = estimate__)) + geom_ribbon(aes(ymin = lower__, ymax = upper__), colour = NA, fill = "#DFBF2B", alpha = .25) + geom_line(colour = "#DFBF2B", linewidth = 2) + scale_y_continuous(limits = c(0,1), breaks = c(0,1)) + labs(x = "Maximum Length (cm)", y = "Probability") + theme_classic() + theme(text = element_text(size = 6), axis.text.x = element_text(angle = 30, hjust = .5, vjust = .75)) ## Trophic Level p07 <- conditional_effects(m01, effects = "Troph", re_formula = NA, spaghetti = FALSE)$Troph %>% ggplot(aes(x = Troph, y = estimate__)) + geom_ribbon(aes(ymin = lower__, ymax = upper__), colour = NA, fill = "#E5A208", alpha = .25) + geom_line(colour = "#E5A208", linewidth = 2) + scale_y_continuous(limits = c(0,1), breaks = c(0,1)) + labs(x = "Trophic Level", y = "Probability") + theme_classic() + theme(text = element_text(size = 6), axis.text.x = element_text(angle = 30, hjust = .5, vjust = .75)) ## Body Shape p08 <- conditional_effects(m01, effects = "BodyShapeI", re_formula = NA, spaghetti = FALSE)$BodyShapeI %>% ggplot(aes(x = reorder(BodyShapeI, -estimate__), y = estimate__)) + geom_pointrange(aes(ymin = lower__, ymax = upper__), colour = "#EA8005", alpha = .25, linewidth = 2) + geom_point(colour = "#EA8005", size = 2) + scale_x_discrete(labels = c("Other", "Eel-like", "Fusiform", "Short/Deep", "Elongated")) + scale_y_continuous(limits = c(0,1), breaks = c(0,1)) + labs(x = "Body Shape", y = "Probability") + theme_classic() + theme(text = element_text(size = 6), axis.text.x = element_text(angle = 30, hjust = .5, vjust = .75)) ## Demersal/Pelagic p09 <- conditional_effects(m01, effects = "DemersPelag", re_formula = NA, spaghetti = FALSE)$DemersPelag %>% ggplot(aes(x = DemersPelag, y = estimate__)) + geom_pointrange(aes(ymin = lower__, ymax = upper__), colour = "#EE5A03", alpha = .25, linewidth = 2) + geom_point(colour = "#EE5A03", size = 2) + scale_y_continuous(limits = c(0,1), breaks = c(0,1)) + scale_x_discrete(labels = c("Benthic", "Benthopelagic", "Pelagic")) + labs(x = "Position", y = "Probability") + theme_classic() + theme(text = element_text(size = 6), axis.text.x = element_text(angle = 30, hjust = .5, vjust = .75)) ## Fertilization p10 <- conditional_effects(m01, effects = "Fertilization", re_formula = NA, spaghetti = FALSE)$Fertilization %>% ggplot(aes(x = Fertilization, y = estimate__)) + geom_pointrange(aes(ymin = lower__, ymax = upper__), colour = "#F11B00", alpha = .25, linewidth = 2) + geom_point(colour = "#F11B00", size = 2) + scale_y_continuous(limits = c(0,1), breaks = c(0,1)) + scale_x_discrete(labels = c("External", "Internal")) + labs(x = "Fertilisation Mode", y = "Probability") + theme_classic() + theme(text = element_text(size = 6), axis.text.x = element_text(angle = 30, hjust = .5, vjust = .75)) ## Plot together pBM <- ggdraw() + draw_plot(plot_grid(p01, p02, p03, p04, p05, p06, p07, p08, p09, p10, ncol = 5, align = "hv"), x = 0, y = 0, width = 1, height = 1) pBM ## clean up rm(p01, p02, p03, p04, p05, p06, p07, p08, p09, p10) gc() ``` ## Counterfactual Now we do our counter-factual predictions Create data for counterfactual prediction ```{r} ## Extract data for counterfactuals cf_final <- complete(d_imp, action = "long") ## back transform the data from the imputation cf_final <- cf_final %>% mutate(sightings = 10^sightings -1, DepthMin = 10^DepthMin -1, DepthMax = 10^DepthMax -1, MaxLengthTL = 10^MaxLengthTL -1, K = 10^K -1, med_depth = 10^med_depth -1, depth_range = 10^depth_range -1, lat_range = 10^lat_range -1 ) ## filter to relevant species only cf_final <- cf_final %>% filter(AS_SC == 1) %>% mutate(across(where(is.character), as.factor)) ## find 95% quantile of sightings - we don't want to set the counterfactual to the most extreme value! s_max <- quantile(cf_final$sightings, 0.95, na.rm = TRUE) ## build dataset for counterfactuals cf_final <- cf_final %>% mutate(sightings_cf = s_max, sightings = sightings_cf) ## predict under counterfactual set.seed(2009) cf_pred <- cf_final %>% add_epred_draws(m01, ndraws = 1000, allow_new_levels = TRUE) %>% ungroup() %>% dplyr::select(c(Phylum:Species, Red_Sea, .epred)) ## mean under counterfactual species_cf <- cf_pred %>% group_by(Class, Order, Family, Genus, Species, Red_Sea) %>% summarise( mean_epred = mean(.epred, na.rm = TRUE), .groups = "drop" ) rm(cf_final) ``` Plot counterfactuals ```{r} ## non-Red Sea species only species_order <- cf_pred %>% filter(Red_Sea == 0) %>% group_by(Class, Order, Family, Genus, Species, Red_Sea) %>% summarise( mean_epred = mean(.epred, na.rm = TRUE), .groups = "drop" ) %>% mutate( Class = factor(Class, levels = c( "Myxini", "Elasmobranchii", "Holocephali", "Teleostei" )) ) %>% arrange(Class, Order, mean_epred) cf_pred <- cf_pred %>% filter(Red_Sea == 0) %>% mutate( Species = factor( Species, levels = species_order$Species ) ) class_labels <- cf_pred %>% distinct(Species, Class) %>% mutate(y = as.numeric(Species)) %>% group_by(Class) %>% summarise( y = mean(y), .groups = "drop" ) class_ranges <- cf_pred %>% distinct(Species, Class) %>% mutate(y = as.numeric(Species)) %>% group_by(Class) %>% summarise( y_min = min(y) - 0.1, y_max = max(y) + 0.1, y_mid = mean(y), .groups = "drop" ) x_outside <- min(cf_pred$.epred) - 0.15 * diff(range(cf_pred$.epred)) x_class <- min(cf_pred$.epred) - 0.1 * diff(range(cf_pred$.epred)) p01 <- ggplot(cf_pred, aes(x = .epred, y = Species)) + geom_density_ridges(rel_min_height = 0.05, colour = "black", scale = 3, bandwidth = 0.05) + geom_text( data = class_labels, aes(x = x_outside, y = y, label = Class), inherit.aes = FALSE, hjust = 1, fontface = "bold", size = 6/.pt ) + ## CLASS RANGE LINES geom_segment( data = class_ranges, aes( x = x_class, xend = x_class, y = y_min, yend = y_max ), inherit.aes = FALSE, linewidth = 0.6 ) + scale_x_continuous( breaks = 0.5, labels = "0.5" ) + geom_vline(xintercept = .5, linetype = 2) + coord_cartesian(clip = "off") + theme_void() + theme( text = element_text(size = 6), legend.position = "none", axis.text.y = element_blank(), axis.ticks.y = element_blank(), plot.margin = margin(5.5, 5.5, 5.5, 45), # LEFT margin space axis.line.x = element_blank(), axis.ticks.x = element_line(), axis.text.x = element_text(size = 6), axis.title.x = element_blank() ) p01 ``` Summarise species and filter high probability (95% CI > .5) and lower probability (50% CI > 0.5) species ```{r} ## summarise CI and demarcate high probability and lower confidence species for inclusion sp_new <- cf_pred %>% group_by(Phylum, Class, Order, Family, Genus, Species) %>% summarise(med = median(.epred), ci2.5 = quantile(.epred, 0.025), ci25 = quantile(.epred, 0.25), ci75 = quantile(.epred, 0.75), ci97.5 = quantile(.epred, 0.975)) %>% mutate(vh_con = ifelse(ci2.5 > 0.5, 1, 0), h_con = ifelse(ci25 > 0.5, 1, 0), l_con = ifelse(med > 0.5, 1, 0)) write.csv(sp_new, "species_predictions.csv", row.names = FALSE) rm(cf_pred) ``` Import image of new species expected ```{r} p02 <- "New_Species/Slide1.png" ``` Plot all together ```{r} pNS <- ggdraw() + draw_plot(pBM, x = .05, y = .7, width = .95, height = .3) + draw_plot(p01, x = 0, y = 0, width = .2, height = .7) + draw_image(p02, x = .2, y = 0, width = .8, height = .7, scale = 1) + draw_plot_label(c("a", "b", "c"), x = c(.01, .01, .21), y = c(.99, .69, .69), size = 8) pNS ggsave("RS_NewSpecies.png", pNS, width = 18, height = 21.5, units = "cm", dpi = 600, bg = "white") ## clean up rm(p01, p02) rm(species_order, class_labels, class_ranges, x_class, x_outside) ``` # Analysis: Trait Space Now we want to examine the trait space that species in the Red Sea inhabit relative to those outside and in the deep sea more generally. ## PCoA Create long form of all the imputations to use for PCoA. Keep key functional traits. Traits should usually cover the below basic areas, listed next to them are traits which are relevant to each: - Food Aquisition: Troph, Feed Pathway, Mouth Position, MaxLengthTL - Mobility: BodyShapeI, MaxLengthTL, lat_range, med_depth, depth_range, DemersPelag - Nutrient Budget: MaxLengthTL, K - Reproduction: Fertilization, MaxLengthTL - Defense Against Predation: Gregariousness, MaxLengthTL We don't use gregariousness values as these available for only 75 of 6707 species We don't use position of mouth values as these available for only 807 of 6707 species Compute the average for each species' traits and use those for modelling: ```{r} ## redraw the dataset, remember that skewed variables here are already log10(var+1) transformed d_PCoA <- complete(d_imp, action = "long") ## keep only traits we want d_PCoA <- d_PCoA %>% select(Species, Troph, MaxLengthTL, FeedingPath, BodyShapeI, lat_range, med_depth, depth_range, DemersPelag, K, Fertilization ) sapply(d_PCoA, class) ## average across traits variables trait_vars <- c("Troph", "MaxLengthTL", "FeedingPath", "BodyShapeI", "lat_range", "med_depth", "depth_range", "DemersPelag", "K", "Fertilization") ## get trait averages d_PCoA <- d_PCoA %>% group_by(Species) %>% summarise( across(all_of(trait_vars), ~ if (is.numeric(.)) mean(.) else names(sort(table(.), decreasing = TRUE))[1] ), .groups = "drop" ) ## convert all characters to factor d_PCoA <- d_PCoA %>% mutate(across(where(is.character), as.factor)) ``` Check variables for skew (just in case) - most are already log10(var+1) transfomed ```{r} hist(d_PCoA$Troph) hist(d_PCoA$MaxLengthTL) hist(d_PCoA$lat_range) hist(d_PCoA$med_depth) hist(d_PCoA$depth_range) hist(d_PCoA$K) colSums(is.na(d_PCoA)) ``` Compute Gower distances ```{r} gower_dist <- daisy( d_PCoA %>% dplyr::select(-Species), metric = "gower" ) ``` Run PCoA ```{r} # Run PCoA once (full eigenvectors) pcoa01 <- ape::pcoa(gower_dist, correction = "cailliez") coords_full <- pcoa01$vectors # all eigenvectors # Prepare vector to store AUCs auc_vals <- numeric(8) # Loop over 1 to 8 axes for(k in 1:8){ message("Calculating AUC using ", k, " axes...") # Select the first k axes coords <- coords_full[, 1:k] # Compute Euclidean distances in reduced space reduced_dist <- dist(coords) # Compute co-ranking matrix Q <- coranking(as.matrix(gower_dist), as.matrix(reduced_dist)) # Compute R_NX curve Rnx <- R_NX(Q) # Compute AUC of R_NX auc_vals[k] <- AUC_ln_K(Rnx) message("Done: AUC = ", round(auc_vals[k], 3)) } # Optional: plot AUC vs number of axes plot(1:8, auc_vals, type="b", xlab="Number of PCoA axes", ylab="AUC", main="AUC vs PCoA axes") ## clean up rm(coords_full, Q, k, reduced_dist, Rnx, coords, gower_dist) gc() ``` Six PCoA axes are required to reach AUC > 0.7 (0.716) Append PCOA scores to the dataframe of species, ensure all regional filter columns are in place ```{r} ## extract PCoA scores scores <- data.frame(pcoa01$vectors[, 1:6]) # give descriptive column names colnames(scores) <- paste0("PCoA", 1:6) ## build species dataframe to append to pcoa_scores <- sp_all %>% left_join(sp_new %>% dplyr::select(Phylum:Species, med)) %>% mutate(RS_pred = ifelse(!is.na(med) & med > 0.5 | Red_Sea == 1, 1, 0)) %>% # Optionally remove the med column if you don’t need it select(-med) %>% dplyr::select(c(Species, IO, AS_SC, RS_pred, Red_Sea)) # Append to species dataframe pcoa_scores <- cbind(pcoa_scores, scores) ## clean rm(scores) ``` Calculate hypervolume of global, IO, AS_SC, and RS spaces and compare to estimate the % truncation of the Red Sea. ```{r} # Compute exact volume for 6D PCoA axes # 'FA' returns the total area and volume hv <- convhulln(pcoa_scores[, 6:11], options = "FA") vol_total <- hv$vol hv1 <- convhulln((pcoa_scores %>% filter(IO == 1))[, 6:11], options = "FA") vol_total1 <- hv1$vol hv2 <- convhulln((pcoa_scores %>% filter(AS_SC == 1))[, 6:11], options = "FA") vol_total2 <- hv2$vol hv3 <- convhulln((pcoa_scores %>% filter(RS_pred == 1))[, 6:11], options = "FA") vol_total3 <- hv3$vol (vol_total1 / vol_total) * 100 ## IO:Global (vol_total2 / vol_total) * 100 ## AS_SC:Global (vol_total3 / vol_total) * 100 ## RS_pred:Global (vol_total3 / vol_total2) * 100 ## RS_pred:AS_SC ``` Correlate variables with PCoA values to see what are the major elements of each axis ```{r} ## traits traits <- d_PCoA[, 2:11] # PCoA axes axes <- pcoa01$vectors[, 1:6] # 1. Initialize empty matrices to store r2 and p-values # Rows = traits, Columns = PCoA Axes axis_names <- paste0("PCoA.", 1:6) trait_names <- colnames(traits) r2_matrix <- matrix(NA, nrow = length(trait_names), ncol = 6, dimnames = list(trait_names, axis_names)) p_matrix <- matrix(NA, nrow = length(trait_names), ncol = 6, dimnames = list(trait_names, axis_names)) # 2. Loop through each axis one by one for (i in 1:6) { # Run envfit on a single column matrix single_axis <- axes[, i, drop = FALSE] fit_single <- envfit(single_axis, traits, permutations = 999) # Extract continuous vectors data if present if (!is.null(fit_single$vectors)) { v_rows <- names(fit_single$vectors$r) r2_matrix[v_rows, i] <- fit_single$vectors$r p_matrix[v_rows, i] <- fit_single$vectors$pvals } # Extract categorical factors data if present if (!is.null(fit_single$factors)) { f_rows <- names(fit_single$factors$r) r2_matrix[f_rows, i] <- fit_single$factors$r p_matrix[f_rows, i] <- fit_single$factors$pvals } } # 3. View the clean R2 breakdown per axis print("--- R2 Values Per Axis ---") print(round(r2_matrix, 4)) print("--- P-Values Per Axis ---") print(round(p_matrix, 4)) ## plot r2 contributions as a heatmap/tile plot library(tidyverse) # Reshape r2_matrix to long format pcoa_r2 <- as.data.frame(r2_matrix) %>% rownames_to_column("Trait") %>% pivot_longer(-Trait, names_to = "Axis", values_to = "R2") %>% mutate(Axis = factor(Axis, levels = paste0("PCoA.", 1:6))) # Plot p01 <- ggplot(pcoa_r2, aes(x = Axis, y = Trait, fill = R2)) + geom_tile(colour = "white", linewidth = 0.5) + scale_y_discrete(labels = c("Body Shape", "Water Column Position", "Depth Range", "Feeding Pathway", "Reproductive Strategy", "K", "Latitudinal Range", "Maximum Length", "Median Depth", "Trophic Level")) + scale_x_discrete(labels = c("PCoA1", "PCoA2", "PCoA3", "PCoA4", "PCoA5", "PCoA6")) + scale_fill_viridis_c(option = "mako", direction = -1, limits = c(0, NA), name = expression(R^2)) + labs(x = NULL, y = NULL) + coord_equal() + theme_minimal(base_size = 12) + theme( axis.text.x = element_text(angle = 45, hjust = 1), panel.grid = element_blank(), text = element_text(size = 6), legend.key.height = unit(0.2, "cm"), legend.key.width = unit(.15, "cm"), legend.position = "bottom", legend.text = element_text(angle = 90, hjust = 1), legend.margin = margin(t = -5, r = 0, b = 0, l = 0, unit = "pt") ) rm(axis_names, trait_names) ``` ## Plot PCoA outputs Plot with convex hulls for subsets of interest ```{r} # Custom panel function — receives the full data + column mappings custom_panel <- function(data, mapping, ...) { xvar <- as_label(mapping$x) yvar <- as_label(mapping$y) hull_all <- data %>% slice(chull(.data[[xvar]], .data[[yvar]])) hull_io <- data %>% filter(IO == 1) %>% slice(chull(.data[[xvar]], .data[[yvar]])) hull_as <- data %>% filter(AS_SC == 1) %>% slice(chull(.data[[xvar]], .data[[yvar]])) hull_rs <- data %>% filter(Red_Sea == 1) %>% slice(chull(.data[[xvar]], .data[[yvar]])) ggplot(data, aes(x = .data[[xvar]], y = .data[[yvar]])) + geom_polygon(data = hull_all, fill = "grey90", color = "grey90") + geom_polygon(data = hull_io, fill = "#899DA4", color = "transparent") + geom_polygon(data = hull_as, fill = "#DC863B", color = "transparent") + geom_polygon(data = hull_rs, fill = "#C93312", color = "black") + geom_point(alpha = 0.1, stroke = 0, size = 0.5, colour = "black") + xlab(xvar) + ylab(yvar) + theme_classic() + theme( panel.border = element_rect(colour = "black", fill = NA, linewidth = 0.4), axis.line = element_blank(), text = element_text(size = 6) ) } p02 <- ggpairs( pcoa_scores, columns = which(names(pcoa_scores) %in% paste0("PCoA", 1:6)), lower = list(continuous = custom_panel), upper = list(continuous = "blank"), diag = list(continuous = "blankDiag")) + coord_cartesian(clip = "on") + theme(strip.background = element_blank(), strip.text = element_blank()) p02 ## extract as a grob (because ggpairs is hard to work with in cowplot) p03 <- grid::grid.grabExpr(print(p02), wrap = TRUE) p03 ``` ## Plot taxonomic trees Plot taxonomic tree with propagation ```{r} # Shared tree tax_tree2 <- as.phylo.formula(~ Phylum/Class/Order/Family/Genus/Species, data = sp_all) tax_tree2 <- compute.brlen(tax_tree2, 1) tax_tree2 <- force.ultrametric(tax_tree2, method = "extend") # Helper: propagate focal status up to ancestors propagate_focal <- function(tree_data, focal_label) { # Mark tips as focal or not tree_data$is_focal <- !is.na(tree_data$edge_region) & tree_data$edge_region == focal_label # Iterate bottom-up: if any child is focal, parent becomes focal # Repeat until no further changes (handles arbitrary tree depth) changed <- TRUE while (changed) { focal_nodes <- tree_data$node[tree_data$is_focal] focal_parents <- unique(tree_data$parent[tree_data$node %in% focal_nodes]) new_focal <- setdiff(focal_parents, focal_nodes) changed <- length(new_focal) > 0 tree_data$is_focal[tree_data$node %in% new_focal] <- TRUE } tree_data$edge_region <- ifelse(tree_data$is_focal, focal_label, "Global") tree_data$is_focal <- NULL tree_data } make_region_tree <- function(tree, sp_data, region_col, region_label, region_hex) { tip_data <- sp_data %>% transmute(label = Species, region = ifelse(.data[[region_col]] == 1, region_label, "Global")) %>% distinct() p <- ggtree(tree, layout = "circular", size = 0) %<+% tip_data p$data <- p$data %>% mutate(edge_region = ifelse(isTip, region, NA)) p$data <- propagate_focal(p$data, region_label) # Full data for both layers; focal edges flagged for top layer d_all <- p$data d_focal <- p$data %>% mutate(edge_region = ifelse(edge_region == region_label, region_label, NA)) p + # Layer 1: full tree in grey (no subsetting, so geometry is intact) geom_tree(data = d_all, colour = "grey90", size = 0.05) + # Layer 2: full data again, but only focal edges get colour (NA = transparent) geom_tree(data = d_focal, aes(colour = edge_region), size = 0.1) + scale_colour_manual( values = setNames(region_hex, region_label), na.value = NA, # non-focal edges in layer 2 are fully transparent guide = "none" ) + theme_tree() + theme(plot.background = element_rect(fill = "transparent", colour = NA), panel.background = element_rect(fill = "transparent", colour = NA), plot.margin = margin(0, 0, 0, 0)) } # Three trees p04 <- make_region_tree(tax_tree2, sp_all, "IO", "IO", "#899DA4") p04 p05 <- make_region_tree(tax_tree2, sp_all, "AS_SC", "AS_SC", "#DC863B") p05 p06 <- make_region_tree(tax_tree2, sp_all, "Red_Sea", "Red Sea", "#C93312") p06 ``` ## Combine plots Plot PcoA and Taxonomic Trees together ```{r} pFS <- ggdraw() + draw_plot(p04, x = 0, y = .64, height = 1/3, width = 1/4) + draw_plot(p05, x = 0, y = .32, height = 1/3, width = 1/4) + draw_plot(p06, x = 0, y = .01, height = 1/3, width = 1/4) + draw_plot(p03, x = 1/4, y = .01, height = 1.15, width = .9) + draw_plot_label(c("PCoA1", "PCoA2", "PCoA3", "PCoA4", "PCoA5"), x = c(.335, .480, .620, .765, .905), y = c(0.01, 0.01, 0.01, 0.01, 0.01), size = 6, hjust = 0, vjust = 0, fontface = "plain") + draw_plot_label(c("PCoA6", "PCoA5", "PCoA4", "PCoA3", "PCoA2"), x = c(.265, .265, .265, .265, .265), y = c(.105, .29, .475, .66, .845), size = 6, hjust = 0, vjust = 0, angle = 90, fontface = "plain") + draw_plot_label(c("Indian Ocean", "Arabian Seas & Somali Current", "Red Sea"), x = c(.125, .125, .125), y = c(.95, .635, .32), colour = c("#899DA4", "#DC863B", "#C93312"), size = 6, hjust = .5, vjust = 0) + draw_plot(p01, x = .69, y = .51, height = .5, width = 1/3) + draw_plot_label(c("a", "b", "c"), x = c(.01, .251, .72), y = c(.99, .99, .99), size = 8) pFS ggsave("RS_Function.png", pFS, width = 16, height = 12, units = "cm", dpi = 600, bg = "white") ```