tutorial.Rmd

---
title: "FISHGLOB Tutorial on Biotic Homogenization"
author: "Malin Pinsky"
date: "2025-12-06"
output: html_document
---

## 1. Introduction
Welcome! In this tutorial, we will explore downloading and using FISHGLOB scientific survey data on fish biodiversity to understand changes in biotic homogenization around the world's oceans. This tutorial is based on the paper by Kitchel et al. 2025 PLOS Climate (https://doi.org/10.1371/journal.pclm.0000659). Our goal today is to re-create the left-hand side of Fig. 2 in Kitchel et al. 2025 for one ocean region.

Many studies report that biological communities are experiencing an increase in spatial similarity through a process termed biotic homogenization. Homogenization can transform overall ecosystem function, lead to instability, and is often considered to be a widespread phenomenon. However, most research to date has focused on terrestrial and freshwater realms, with little understanding for marine ecosystems. The FISHGLOB dataset provides a first window into biotic homogenization around the world for marine continental shelf ecosystems.

To proceed through this tutorial, read the text above each grey chunk of code. To run the code in a chunk, either click the green arrow in the upper right corner of a chunk, or select some lines of code and choose "Run Selected Lines" from the Code menu.

This is some setup code. Please run it now. If you don't have the "vegan" or "ggplot2" packages, it will install them for you.
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)

# install packages that are not on your machine yet
if(!require("vegan")) install.packages("vegan")
if(!require("ggplot2")) install.packages("ggplot2")

# load packages once installed
library(vegan) # for calculating dissimilarity in species compositions
library(ggplot2) # for plotting
```

## 2. Download data
Our first step is to download the FISHGLOB data. The data is stored in a public website (technically a "GitHub repository") here: https://github.com/fishglob/FishGlob_data along with R code to assemble the data. 

The cleaned data are stored on the website under outputs/Cleaned_data: https://github.com/fishglob/FishGlob_data/tree/main/outputs/Cleaned_data. Each region has an abbreviation (e.g., NEUS for the Northeast United States).

To download the data, all we need to do is point R to the address of the file. Let's try the Northeast US first. Run the following code to load the data into R:
```{r load data}
options(timeout = 240) # set a longer download timeout in case you are on a slow connection
load(url("https://github.com/fishglob/FishGlob_data/raw/refs/heads/flag_hex/outputs/Cleaned_data/NEUS_std_clean.RData")) # for the Northeast US
#load(url("https://github.com/fishglob/FishGlob_data/raw/refs/heads/flag_hex/outputs/Cleaned_data/GMEX_std_clean.RData")) # another option, for the Gulf of Mexico.
```

## 3. View data
The survey data are stored as a data table (like an Excel spreadsheet) with many columns in an object called "data." For the description of all the columns, see the descriptions at https://github.com/fishglob/FishGlob_data/tree/main/standard_formats.

Each row in this spreadsheet is an observation of a species at a particular location and time during a haul of a bottom trawl survey.

This code will provide a summary of the dataset columns and a small look at what is in each one. The "survey" column, for example, provides the survey name. The "source" column provides the name of the institution that performed the survey. The "year", "month", and "day" columns provide the date of a particular survey tow. The wgt_cpue and the num_cpue provide observations of biomass (wgt_cpue) and abundance (num_cpue) standardized based on the observation effort.
```{r view}
str(data)
```

For a description of all the columns, look at the "readme" object that is provided with each FISHGLOB survey file. This was loaded with the data.
```{r metadata}
View(readme)
```

We can also plot the locations of the survey tows. A transparent color is used to show which locations were surveyed more often.
```{r map}
plot(data$longitude, data$latitude, col = '#00000001', xlab = "Longitude", ylab = "Latitude")

# code to plot what will be trimmed out (black) or kept (red) in the spatial standardization step later.
# ggplot(data, aes(longitude, latitude)) + 
#     geom_point(color = is.na(data$flag_trimming_hex8_2)+1, size = 0.2) +
#     facet_wrap(vars(year))

```

## 4. Clean up
We need to clean up the dataset to prepare for the analysis of biotic homogenization. Depending on the analysis one is doing, these steps will be different.

First, we will remove the observations not resolved to species or subspecies. The "rank" column provides a handy indication of the taxonomic rank of each entry (e.g., family, genus, species, etc.).
```{r spptrim}
datatrim <- data[data$rank %in% c("Species", "Subspecies"),]
```

Second, we will trim every year to a consistent spatial footprint (see Fig. 3 in Maureaud et al. 2024 https://doi.org/10.1038/s41597-023-02866-w for a visual representation of this step). The "flag_trimming" columns provide some pre-calculated options. We will use "flag_trimming_8_0", which is the finer spatial grain option and allows 0% missing years from each spatial grid cell. TRUE indicates a row to remove, NA indicates a row to keep.
```{r spatialtrim}
datatrim.spatial <- datatrim[is.na(datatrim$flag_trimming_hex8_0),]
```

Third, delete all rows with missing abundance and biomass data.
```{r NA}
datatrim.noNA <- datatrim.spatial[!(is.na(datatrim.spatial$num_cpue) & is.na(datatrim.spatial$wgt_cpue) & is.na(datatrim.spatial$wgt) & is.na(datatrim.spatial$num)),]
```

Let's plot the locations of the survey tows again, for a comparison against the full dataset we mapped earlier. We can see how some poorly surveyed parts were trimmed off.
```{r map}
plot(datatrim.noNA$longitude, datatrim.noNA$latitude, col = '#00000001', xlab = "Longitude", ylab = "Latitude")

# plot by year
# ggplot(datatrim.noNA, aes(longitude, latitude)) + 
#     geom_point(size = 0.2) +
#     facet_wrap(vars(year))


```

## 5. Calculate homogenization in each year
We use the Jaccard dissimilarity metric among locations in each year to understand the degree of homogenization each year. We make this calculation each year and then can plot the changes through time. This calculation takes a little while to loop through each year of the survey (a few seconds).
```{r jaccard}
#list the years of data we have
years <- unique(datatrim.noNA$year)

# Create an object to hold the dissimilarites
jaccard_index <- data.frame("year" = integer(length(years)),
                            "annual_dissimilarity_value" = numeric(length(years)))


# create a column for whether a species is present (TRUE) or absent (FALSE).
# a species is present if abundance (num and num_cpue) or biomass (wgt and wgt_cpue) is greater than zero
datatrim.noNA$Present <- ifelse(datatrim.noNA$num_cpue>0 | datatrim.noNA$wgt_cpue>0 | datatrim.noNA$wgt>0 | datatrim.noNA$num>0, TRUE, FALSE)

for (j in 1:length(years)) {
    # trim to data for this year
    colstokeep <- c("haul_id", "Present", "accepted_name")
    reduced_year <- unique(datatrim.noNA[datatrim.noNA$year == years[j], colstokeep])
    
    # convert to a data.frame so that reshape() will work properly
    reduced_year <- as.data.frame(reduced_year)
    
    #Only include rows where species are present
    reduced_year_occurence <- reduced_year[reduced_year$Present==TRUE, ]
    
    #Convert species name to a number for easier reshaping
    reduced_year_occurence$accepted_name <- as.numeric(as.factor(reduced_year_occurence$accepted_name))

    # Make a community matrix (rows are years, columns are species)
    reduced_year_wide <- reshape(reduced_year_occurence, direction = "wide", idvar = "haul_id", timevar = "accepted_name") #convert to wide format
    ncols <- ncol(reduced_year_wide)
    communitymatrix <- reduced_year_wide[, 2:ncols] #extract the community matrix (drop the haul_id column)
    communitymatrix[is.na(communitymatrix)] <- 0 # fill 0s for species/location combinations not observed (NA)
    
    #Calculate Jaccard dissimilarity values
    Jaccard_index_binary <- vegdist(communitymatrix,  method = "jaccard", binary = TRUE)

    #take the average and save the value from this year
    jaccard_index$year[j] <- years[j]
    jaccard_index$annual_dissimilarity_value[j] <- mean(Jaccard_index_binary, na.rm=TRUE)
}
```


## 6. Plot
Time to visualize the results! We plot year on the x-axis and dissimilarity on the y-axis (down is more homogenous, up is more dissimilar). We also plot a curve through the datapoints.
```{r plot}

ggplot(jaccard_index, aes(year, annual_dissimilarity_value)) +
    geom_smooth() +
    geom_point(alpha = 0.4, size = 1, shape = 21) +
    theme_classic() +
    xlab("Year") +
    ylab("β-diversity\n(Jaccard dissimilarity)") +
    theme(legend.position = "null", axis.text = element_text(size = 15), axis.title = element_text(size = 15))

```

## 7. Questions

1. For this region, how is the degree of homogenization changing? Is it clearly increasing or decreasing?

2. If you know the marine community in this region, can you think of reasons that the degree of homogenization might be changing?

3. Go back to Step 2 (Download data) and pick a different region to download and visualize. For a map of the options, check out Fig. 2 of Maureaud et al. 2024 (https://doi.org/10.1038/s41597-023-02866-w). You can see the list of files for different regions in the outputs/Cleaned_data folder of the Github repository (https://github.com/fishglob/FishGlob_data/tree/main/outputs/Cleaned_data). Make sure to download the _std_clean.RData version so that the spatial standardization flags are included.