Association of invasive rodents on native species abundance, richness and diversity
Msc thesis- Conservation Management of African Ecosystems
Maya Collins
23-08-2024
Loading necessary packages
# Load necessary packages for data manipulation and visualization
library(lubridate) # For date-time manipulation##
## Attaching package: 'lubridate'## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, unionlibrary(dplyr) # For data manipulation##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(sf) # For spatial data handling## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUElibrary(tidyr) # For data tidying
library(raster) # For raster data handling## Loading required package: sp##
## Attaching package: 'raster'## The following object is masked from 'package:dplyr':
##
## selectlibrary(tidyr) # For data tidying
library(ggplot2) # For data visualisation
library(vegan) # For ecological data analysis## Loading required package: permute## Loading required package: lattice## This is vegan 2.6-6.1library(ggeffects) # For plotting regression effects
library(MASS) # For negative binomial regression##
## Attaching package: 'MASS'## The following objects are masked from 'package:raster':
##
## area, select## The following object is masked from 'package:dplyr':
##
## selectlibrary(DHARMa) # For checking model assumptions## This is DHARMa 0.4.6. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa')library(terra) # For advanced raster data handling## terra 1.7.78##
## Attaching package: 'terra'## The following object is masked from 'package:MASS':
##
## area## The following object is masked from 'package:tidyr':
##
## extractlibrary(gridExtra) # For arranging multiple plots##
## Attaching package: 'gridExtra'## The following object is masked from 'package:dplyr':
##
## combinelibrary(car) # For checking multicollinearity## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recodelibrary(broom) # For tidying model output
library(kableExtra) # For creating HTML tables##
## Attaching package: 'kableExtra'## The following object is masked from 'package:dplyr':
##
## group_rowsData Cleaning, Organising and Manipulation
Reading rodent data
This data includes species information, geographical coordinates, and the year of capture.
# Read in rodent data from CSV file
data_GBIF<-read.csv("/Users/mayacollins/Desktop/Msc Research project/GBIF_2000-2017.csv", header = T)Cleaning and Filtering Data
# Count the number of empty strings in species names
empty_species_count <- sum(data_GBIF$species == "")
cat("Number of empty species names:", empty_species_count, "\n")## Number of empty species names: 744# Filter out rows with empty species names
data_GBIF_clean <- data_GBIF %>%
filter(species != "")Creating Trap events
# Create unique trap event identifier for rodents based on year, longitude, and latitude
data_GBIF$trap_event_2<-paste(data_GBIF$year,data_GBIF$Longitude, data_GBIF$Latitude, sep = "_")
# Calculate the total number of unique trap events
num_trap_events <- length(levels(as.factor(data_GBIF$trap_event_2)))
cat("Total number of trap events:", num_trap_events, "\n") ## Total number of trap events: 390# Calculate average number of trapped rodents per trap event
avg_rodents_per_event <- nrow(data_GBIF) / num_trap_events
cat("Average number of trapped rodents per trap event:", avg_rodents_per_event, "\n")## Average number of trapped rodents per trap event: 30.7641Categorising invasive species
# Define invasive species
invasive_species <- c('Rattus rattus', 'Rattus norvegicus', 'Mus musculus')
# Add a column indicating species status (invasive or non-invasive)
data_GBIF <- data_GBIF %>%
mutate(
species_status = ifelse(species %in% invasive_species, 'invasive', 'non-invasive')
)Importing and visualising protected areas
# Importing protected areas shapefile of Uganda, Kenya and Tanzania
PA <- st_read("/Users/mayacollins/Desktop/Msc Research project/shapefiles/EA.shp")## Reading layer `EA' from data source
## `/Users/mayacollins/Desktop/Msc Research project/shapefiles/EA.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 625 features and 32 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 29.5836 ymin: -11.28136 xmax: 41.5642 ymax: 4.260569
## Geodetic CRS: WGS 84# Plot the protected areas to visually inspect their geographical distribution
ggplot() +
geom_sf(data = PA, size = 1.5, color = "black", fill = "cyan1") +
ggtitle("East Africa PAs") +
coord_sf()
Creating a New Dataset with Trap Events
# Summarize the data by trap event to calculate total individuals, species richness, and invasive/native abundance
trap_summ <- data_GBIF %>%
group_by(trap_event_2, year) %>%
summarise(
total_individuals = n(),
species_list = list(unique(species)),
species_list_native = list(unique(species[species_status == 'non-invasive'])),
invasive_abundance = sum(species_status == 'invasive'),
native_abundance = sum(species_status == 'non-invasive'),
native_species_richness = n_distinct(species[species_status == 'non-invasive']),
latitude = first(Lattitude),
longitude = first(Longitude),
relative_abundance_invasive = sum(species_status =='invasive') / n(),
relative_abundance_native = sum(species_status =='non-invasive') / n() )## `summarise()` has grouped output by 'trap_event_2'. You can override using the
## `.groups` argument.Adding Protected Area Information to Trap Events
# Convert trap summary data to a spatial object for spatial analysis
trap_summ_sf <- st_as_sf(trap_summ, coords = c("longitude", "latitude"), crs = st_crs(PA))
# Check if points are inside protected areas
trap_summ_sf$inside_PA <- ifelse(lengths(st_intersects(trap_summ_sf, PA)) > 0, "inside", "outside")
# Convert back to data frame
trap_summ <- as.data.frame(trap_summ_sf)
# Add latitude and longitude back to data frame
trap_summ$latitude <- st_coordinates(trap_summ_sf)[, 2]
trap_summ$longitude <- st_coordinates(trap_summ_sf)[, 1]Calculating the Percentage of Protected Areas with Invasive Species
# Filter traps with invasive species
traps_with_invasive <- trap_summ_sf %>% filter(invasive_abundance > 0)
# Check if each trap with invasive species is within any protected area
inside_PA_invasive <- st_within(traps_with_invasive, PA, sparse = FALSE)
# Get the indices of protected areas that have invasive species
pa_with_invasive_indices <- which(rowSums(inside_PA_invasive) > 0)
# Extract the corresponding protected areas
pa_with_invasive <- PA[pa_with_invasive_indices, ]
# Calculate the percentage of protected areas with invasive species
percentage_with_invasive <- (nrow(pa_with_invasive) / nrow(PA)) * 100
# Print the result
percentage_with_invasive## [1] 11.2Determining if each protected area has invasive species
# Determine if invasive species are present based on invasive abundance at each trap event
trap_summ <- trap_summ %>%
mutate(PA_invasive_presence = ifelse(invasive_abundance > 0, "with_invasive", "without_invasive"))Visualising trap events in relation to protected areas
# Visualize the trap events in relation to the protected areas to explore their spatial distribution
ggplot() +
geom_sf(data = PA, size = 1.5, color = "black", fill = "forestgreen") +
geom_point(data = trap_summ, aes(x = longitude, y = latitude, color = inside_PA), size = 1) +
coord_sf() +
theme_minimal() +
scale_color_manual(values = c("inside" = "hotpink", "outside" = "blue"), name = "Trap Location")
Calculate the proportions of traps with invasive
# Filter the dataset to find trap events with invasive species
trap_with_invasive <- trap_summ %>%
filter(invasive_abundance > 0)
# Calculate the number of trap events with invasive species
num_trap_with_invasive <- nrow(trap_with_invasive)
# Calculate the total number of trap events
total_trap_events <- nrow(trap_summ)
# Calculate the percentage of trap events with invasive species
percentage_with_invasive <- (num_trap_with_invasive / total_trap_events) * 100
# Print the results
cat("Number of trap events with invasive species:", num_trap_with_invasive, "\n")## Number of trap events with invasive species: 105cat("Total number of trap events:", total_trap_events, "\n")## Total number of trap events: 390cat("Percentage of trap events with invasive species:", percentage_with_invasive, "%", "\n")## Percentage of trap events with invasive species: 26.92308 %Identify unique protected areas with trap events
# Find unique protected areas that have any trap events
pa_with_trap_events <- st_intersects(trap_summ_sf, PA, sparse = FALSE)
# Get unique protected areas that intersect with any trap event
pa_with_trap_events_indices <- unique(which(rowSums(pa_with_trap_events) > 0))
# Total number of protected areas with trap events
num_pa_with_trap_events <- length(pa_with_trap_events_indices)Identify protected areas with invasive species and calculate perecentages
# Find unique protected areas that have invasive species
pa_with_invasive_trap_events <- st_intersects(traps_with_invasive, PA, sparse = FALSE)
# Get unique protected areas that intersect with any invasive species trap event
pa_with_invasive_trap_events_indices <- unique(which(rowSums(pa_with_invasive_trap_events) > 0))
# Number of protected areas with invasive species
num_pa_with_invasive <- length(pa_with_invasive_trap_events_indices)
# Total number of protected areas
total_pa <- nrow(PA)
# Percentage of protected areas that have trap events in them
percentage_pa_with_trap_events <- (num_pa_with_trap_events / total_pa) * 100
# Percentage of protected areas with trap events that have invasive species in them
percentage_pa_with_invasive <- (num_pa_with_invasive / num_pa_with_trap_events) * 100
# Print the results
cat("Total number of protected areas:", total_pa, "\n")## Total number of protected areas: 625cat("Number of protected areas with trap events:", num_pa_with_trap_events, "\n")## Number of protected areas with trap events: 242cat("Percentage of protected areas with trap events:", percentage_pa_with_trap_events, "%\n")## Percentage of protected areas with trap events: 38.72 %cat("Number of protected areas with invasive species:", num_pa_with_invasive, "\n")## Number of protected areas with invasive species: 70cat("Percentage of protected areas with invasive species:", percentage_pa_with_invasive, "%\n")## Percentage of protected areas with invasive species: 28.92562 %Calculate Shannon Diversity Index for Each Trap Event
# Filter data to include only native species
native_species_data <- data_GBIF %>%
filter(species_status == 'non-invasive')
# Remove rows with empty or NA species names
native_species_data <- native_species_data %>%
filter(!is.na(species) & species != "")
# Create a species abundance matrix (rows are trap events, columns are species)
species_abundance_matrix <- native_species_data %>%
group_by(trap_event_2, species) %>%
summarise(abundance = n(), .groups = 'drop') %>%
pivot_wider(names_from = species, values_from = abundance, values_fill = list(abundance = 0))
# Set trap_event_2 as row names
species_abundance_matrix <- as.data.frame(species_abundance_matrix)
rownames(species_abundance_matrix) <- species_abundance_matrix$trap_event_2
species_abundance_matrix <- species_abundance_matrix[ , -1]
# Calculate Shannon diversity index for each trap event
shannon_diversity <- diversity(species_abundance_matrix, index = "shannon")
# Convert the result to a data frame
shannon_diversity_df <- data.frame(trap_event_2 = rownames(species_abundance_matrix), shannon_diversity = shannon_diversity)Merge with Trap Summary Data
# Merge Shannon diversity index with trap summary data
trap_summ <- trap_summ %>%
left_join(shannon_diversity_df, by = "trap_event_2")Converting columns to factors
# Convert relevant columns to factor type
trap_summ$inside_PA <- factor(trap_summ$inside_PA)
trap_summ$trap_event_2 <- factor(trap_summ$trap_event_2)
trap_summ$PA_invasive_presence <- factor(trap_summ$PA_invasive_presence)
glimpse(trap_summ) # Verify structure of the data## Rows: 390
## Columns: 16
## $ trap_event_2 <fct> 2000_32.8167_, 2000_33.9231_, 2000_35.1667…
## $ year <int> 2000, 2000, 2000, 2000, 2000, 2000, 2000, …
## $ total_individuals <int> 10, 20, 3, 11, 6, 2, 8, 3, 5, 4, 9, 5, 5, …
## $ species_list <list> <"Mus tenellus", "Aethomys chrysophilus",…
## $ species_list_native <list> <"Mus tenellus", "Aethomys chrysophilus",…
## $ invasive_abundance <int> 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ native_abundance <int> 10, 20, 3, 10, 6, 2, 8, 3, 5, 4, 9, 5, 5, …
## $ native_species_richness <int> 3, 3, 1, 4, 4, 2, 6, 3, 3, 4, 4, 1, 3, 1, …
## $ relative_abundance_invasive <dbl> 0.00000000, 0.00000000, 0.00000000, 0.0909…
## $ relative_abundance_native <dbl> 1.0000000, 1.0000000, 1.0000000, 0.9090909…
## $ geometry <POINT [°]> POINT (32.8167 -5.4167), POINT (33.9…
## $ inside_PA <fct> outside, inside, outside, outside, outside…
## $ latitude <dbl> -5.4167, -9.0636, -7.8167, -7.7997, -8.066…
## $ longitude <dbl> 32.8167, 33.9231, 35.1667, 35.7658, 35.900…
## $ PA_invasive_presence <fct> without_invasive, without_invasive, withou…
## $ shannon_diversity <dbl> 1.0889000, 0.7305881, 0.0000000, 1.0889000…Loading MODIS NDVI data from 2000-2017 (1km resolution)
# This function loads NDVI files for each month of a given year and calculates the yearly average NDVI.
# Parameters:
# year: A string representing the year for which NDVI data should be loaded (e.g., "2000").
# Returns:
# ndvi_avg: A raster object representing the average NDVI for the specified year.
load_ndvi_data <- function(year) {
ndvi_path <- paste0("/Users/mayacollins/Desktop/Msc Research project/MODIS_NDVI_Maya/MODIS_NDVI_", year, "/") # # Construct the file path to the directory containing NDVI files for the given year
files <- list.files(ndvi_path) # List all NDVI files in the directory for that year
ndvi_list <- lapply(files, function(file) {
rast(paste0(ndvi_path, file)) # Load each NDVI file into a list as a raster object
})
ndvi_avg <- do.call(mean, c(ndvi_list, list(na.rm = TRUE))) # Calculate the average NDVI across all months, ignoring any NA values
names(ndvi_avg) <- paste0("NDVI_average_", year) # Rename the resulting raster layer to indicate the year
return(ndvi_avg) # Return the average NDVI raster
}
# Calculating the average NDVI for each year from 2000 to 2017
NDVI_avg_2000 <- load_ndvi_data("2000")
NDVI_avg_2001 <- load_ndvi_data("2001")
NDVI_avg_2002 <- load_ndvi_data("2002")
NDVI_avg_2003 <- load_ndvi_data("2003")
NDVI_avg_2004 <- load_ndvi_data("2004")
NDVI_avg_2005 <- load_ndvi_data("2005")
NDVI_avg_2006 <- load_ndvi_data("2006")
NDVI_avg_2007 <- load_ndvi_data("2007")
NDVI_avg_2008 <- load_ndvi_data("2008")
NDVI_avg_2009 <- load_ndvi_data("2009")
NDVI_avg_2010 <- load_ndvi_data("2010")
NDVI_avg_2011 <- load_ndvi_data("2011")
NDVI_avg_2012 <- load_ndvi_data("2012")
NDVI_avg_2013 <- load_ndvi_data("2013")
NDVI_avg_2014 <- load_ndvi_data("2014")
NDVI_avg_2015 <- load_ndvi_data("2015")
NDVI_avg_2016 <- load_ndvi_data("2016")
NDVI_avg_2017 <- load_ndvi_data("2017")Plot and animate
all_NDVI<-rast(list(NDVI_avg_2000,
NDVI_avg_2001,
NDVI_avg_2002,
NDVI_avg_2003,
NDVI_avg_2004,
NDVI_avg_2005,
NDVI_avg_2006,
NDVI_avg_2007,
NDVI_avg_2008,
NDVI_avg_2009,
NDVI_avg_2010,
NDVI_avg_2011,
NDVI_avg_2012,
NDVI_avg_2013,
NDVI_avg_2014,
NDVI_avg_2015,
NDVI_avg_2016,
NDVI_avg_2017))
plot(all_NDVI)
Extracting NDVI and for every year
buff<-st_buffer(trap_summ_sf, dist=5000)
trap_summ$NDVI<-NA
temp_NDVI2000<-extract(NDVI_avg_2000, buff, fun=mean, na.rm=T)
temp_NDVI2001<-extract(NDVI_avg_2001, buff, fun=mean, na.rm=T)
temp_NDVI2002<-extract(NDVI_avg_2002, buff, fun=mean, na.rm=T)
temp_NDVI2003<-extract(NDVI_avg_2003, buff, fun=mean, na.rm=T)
temp_NDVI2004<-extract(NDVI_avg_2004, buff, fun=mean, na.rm=T)
temp_NDVI2005<-extract(NDVI_avg_2005, buff, fun=mean, na.rm=T)
temp_NDVI2006<-extract(NDVI_avg_2006, buff, fun=mean, na.rm=T)
temp_NDVI2007<-extract(NDVI_avg_2007, buff, fun=mean, na.rm=T)
temp_NDVI2008<-extract(NDVI_avg_2008, buff, fun=mean, na.rm=T)
temp_NDVI2009<-extract(NDVI_avg_2009, buff, fun=mean, na.rm=T)
temp_NDVI2010<-extract(NDVI_avg_2010, buff, fun=mean, na.rm=T)
temp_NDVI2011<-extract(NDVI_avg_2011, buff, fun=mean, na.rm=T)
temp_NDVI2012<-extract(NDVI_avg_2012, buff, fun=mean, na.rm=T)
temp_NDVI2013<-extract(NDVI_avg_2013, buff, fun=mean, na.rm=T)
temp_NDVI2014<-extract(NDVI_avg_2014, buff, fun=mean, na.rm=T)
temp_NDVI2015<-extract(NDVI_avg_2015, buff, fun=mean, na.rm=T)
temp_NDVI2016<-extract(NDVI_avg_2016, buff, fun=mean, na.rm=T)
temp_NDVI2017<-extract(NDVI_avg_2017, buff, fun=mean, na.rm=T)Adding NDVI in the trap event dataset
trap_summ$NDVI[which(trap_summ$year==2000)]<-temp_NDVI2000[which(trap_summ$year==2000),2]
trap_summ$NDVI[which(trap_summ$year==2001)]<-temp_NDVI2001[which(trap_summ$year==2001),2]
trap_summ$NDVI[which(trap_summ$year==2002)]<-temp_NDVI2002[which(trap_summ$year==2002),2]
trap_summ$NDVI[which(trap_summ$year==2003)]<-temp_NDVI2003[which(trap_summ$year==2003),2]
trap_summ$NDVI[which(trap_summ$year==2004)]<-temp_NDVI2004[which(trap_summ$year==2004),2]
trap_summ$NDVI[which(trap_summ$year==2005)]<-temp_NDVI2005[which(trap_summ$year==2005),2]
trap_summ$NDVI[which(trap_summ$year==2006)]<-temp_NDVI2006[which(trap_summ$year==2006),2]
trap_summ$NDVI[which(trap_summ$year==2007)]<-temp_NDVI2007[which(trap_summ$year==2007),2]
trap_summ$NDVI[which(trap_summ$year==2008)]<-temp_NDVI2008[which(trap_summ$year==2008),2]
trap_summ$NDVI[which(trap_summ$year==2009)]<-temp_NDVI2009[which(trap_summ$year==2009),2]
trap_summ$NDVI[which(trap_summ$year==2010)]<-temp_NDVI2010[which(trap_summ$year==2010),2]
trap_summ$NDVI[which(trap_summ$year==2011)]<-temp_NDVI2011[which(trap_summ$year==2011),2]
trap_summ$NDVI[which(trap_summ$year==2012)]<-temp_NDVI2012[which(trap_summ$year==2012),2]
trap_summ$NDVI[which(trap_summ$year==2013)]<-temp_NDVI2013[which(trap_summ$year==2013),2]
trap_summ$NDVI[which(trap_summ$year==2014)]<-temp_NDVI2014[which(trap_summ$year==2014),2]
trap_summ$NDVI[which(trap_summ$year==2015)]<-temp_NDVI2015[which(trap_summ$year==2015),2]
trap_summ$NDVI[which(trap_summ$year==2016)]<-temp_NDVI2016[which(trap_summ$year==2016),2]
trap_summ$NDVI[which(trap_summ$year==2017)]<-temp_NDVI2017[which(trap_summ$year==2017),2]
rm(list = ls(pattern = "temp_NDVI")) #removing temporary NDVI filesModeling and Model Prediction Plots
Analying Invasive and Native Species Associations
Native Presence
Objective: Investigate the relationship between invasive rodent abundance species and native rodent abundance, incorporating NDVI (Normalized Difference Vegetation Index) and the location of trap events relative to protected areas (inside/outside).
Splitting the data for model analysis. The analysis will be using two separate models to address the high number of zeros in the data set. The first model is a binary logistic regression to analyse the presence/absence of native rodents, and the second model is a negative binomial regression to analyse the abundance of native rodents where they are present.
# Convert native_abundance to binary presence/absence
trap_summ$native_presence <- ifelse(trap_summ$native_abundance > 0, 1, 0)
# Fit logistic regression model
presence_model <- glm(native_presence ~ invasive_abundance + NDVI + inside_PA, family = binomial, data = trap_summ)
summary(presence_model) # summarising model output##
## Call:
## glm(formula = native_presence ~ invasive_abundance + NDVI + inside_PA,
## family = binomial, data = trap_summ)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.02249 0.92405 -1.107 0.268497
## invasive_abundance -0.05969 0.01656 -3.605 0.000312 ***
## NDVI 5.66230 1.62458 3.485 0.000491 ***
## inside_PAoutside 1.13034 0.48229 2.344 0.019095 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 252.72 on 385 degrees of freedom
## Residual deviance: 200.71 on 382 degrees of freedom
## (4 observations deleted due to missingness)
## AIC: 208.71
##
## Number of Fisher Scoring iterations: 6simms_pr <- simulateResiduals(fittedModel = presence_model) # generates simulated residuals for the logistic regression model presence_model using the simulateResiduals function
plotQQunif(simms_pr) # creates a quantile-quantile (Q-Q) plot of the simulated residuals using the plotQQunif function
vif(presence_model) #checking for multicolleniarity among predictors## invasive_abundance NDVI inside_PA
## 1.068561 1.047840 1.039720tidy(presence_model, conf.int = TRUE) # extracting confidence intervals## # A tibble: 4 × 7
## term estimate std.error statistic p.value conf.low conf.high
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) -1.02 0.924 -1.11 0.268 -2.85 0.800
## 2 invasive_abundance -0.0597 0.0166 -3.61 0.000312 -0.0946 -0.0296
## 3 NDVI 5.66 1.62 3.49 0.000491 2.56 8.98
## 4 inside_PAoutside 1.13 0.482 2.34 0.0191 0.244 2.17Plotting from Native Presence and Abundance Predictions
gg_ndvi <- ggpredict(presence_model, terms = c("NDVI[all]")) # generates predicted values from the logistic regression model presence_model for the term "NDVI" using ggpredict
pres1 <- ggplot(gg_ndvi, aes(x = x, y = predicted)) + # plots data from ggpredict
geom_line() + # adds line to the plot representing the predicted probabilities
geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = 0.2, colour = NA) + #Adds a shaded area (ribbon) around the line to represent the 95% confidence intervals.
theme_bw() + # Applies a clean, white background theme to the plot.
labs(x = "NDVI", y = "Probability of Native Presence", fill = "Confidence Intervals") # Labels the x-axis as "NDVI", the y-axis as "Probability of Native Presence", and the fill legend as "Confidence Intervals".
gg_invasive <- ggpredict(presence_model, terms = c("invasive_abundance[all]")) # generates predicted values from the logistic regression model presence_model for the term "invasive_abundance" using ggpredict
pres2 <- ggplot(gg_invasive, aes(x = x, y = predicted)) +
geom_line() +
geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = 0.2, colour = NA) +
theme_bw() +
labs(x = "Invasive Abundance", y = "Probability of Native Presence", fill = "Confidence Intervals")
grid.arrange(pres1, pres2, ncol=2) # arranges the two plots (pres1 and pres2) side by side in a single output grid, with 2 columns (ncol=2).
Native Abundance
# Subset the data to include only rows where native rodents are present
non_zero_data <- subset(trap_summ, native_abundance > 0)
# Fit negative binomial model on the subset
abundance_model <- glm.nb(native_abundance ~ invasive_abundance + NDVI +inside_PA, data = non_zero_data)
summary(abundance_model)##
## Call:
## glm.nb(formula = native_abundance ~ invasive_abundance + NDVI +
## inside_PA, data = non_zero_data, init.theta = 0.5383635065,
## link = log)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.88419 0.38897 4.844 1.27e-06 ***
## invasive_abundance -0.01844 0.01113 -1.657 0.09762 .
## NDVI 2.72831 0.59156 4.612 3.99e-06 ***
## inside_PAoutside -0.43182 0.15279 -2.826 0.00471 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.5384) family taken to be 1)
##
## Null deviance: 443.13 on 346 degrees of freedom
## Residual deviance: 412.34 on 343 degrees of freedom
## (3 observations deleted due to missingness)
## AIC: 2951.7
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.5384
## Std. Err.: 0.0357
##
## 2 x log-likelihood: -2941.6970simms_ab <- simulateResiduals(fittedModel = abundance_model)
plotQQunif(simms_ab)
vif(abundance_model)## invasive_abundance NDVI inside_PA
## 1.059707 1.056270 1.026594tidy(abundance_model, conf.int = TRUE)## # A tibble: 4 × 7
## term estimate std.error statistic p.value conf.low conf.high
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 1.88 0.389 4.84 0.00000127 1.21 2.59
## 2 invasive_abundance -0.0184 0.0111 -1.66 0.0976 -0.0373 0.00657
## 3 NDVI 2.73 0.592 4.61 0.00000399 1.62 3.82
## 4 inside_PAoutside -0.432 0.153 -2.83 0.00471 -0.732 -0.127Plotting from Native Abundance Predictions
gg_ndvi <- ggeffect(abundance_model, terms = c("NDVI [all]", "inside_PA")) # generates predicted effects of NDVI on native rodent abundance from the abundance_model, distinguishing between inside and outside protected areas.
# Create NDVI plot with inside_PA using ggeffects
abun_ndvi <- ggplot(gg_ndvi, aes(x = x, y = predicted, colour = group)) +
geom_line() +
geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = group), alpha = 0.2, colour = NA) +
theme_minimal() +
labs(x = "NDVI",y = "Native Abundance",colour = "Inside PA",fill = "Confidence Intervals")
gg_invasive_abun <- ggeffect(abundance_model, terms = c("invasive_abundance [all]", "inside_PA")) # generates predicted effects of invasive rodent abundance on native rodent abundance, also distinguishing between inside and outside protected areas.
# Create Invasive Abundance plot with inside_PA using ggeffects
abun_invasive <- ggplot(gg_invasive_abun, aes(x = x, y = predicted, colour = group)) +
geom_line() +
geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = group), alpha = 0.2, colour = NA) +
theme_minimal() +
labs(x = "Invasive Abundance",y = "Native Abundance",colour = "Inside PA", fill = "Confidence Intervals")
grid.arrange(abun_ndvi, abun_invasive, nrow = 2)
Investigate Invasive Invasive Rodent Presence/Abundance Diversity and Species Richness
Native Richness
Objective: Examine how the abundance of invasive rodents affects the diversity and species richness of native rodent species abundance, while considering NDVI and the trap event location (inside vs. outside protected areas).
# Negative binomial model with offset to analyse the relationship
glm_richness <- glm.nb(native_species_richness ~ PA_invasive_presence + inside_PA + NDVI, data = trap_summ)
summary(glm_richness)##
## Call:
## glm.nb(formula = native_species_richness ~ PA_invasive_presence +
## inside_PA + NDVI, data = trap_summ, init.theta = 3.144535995,
## link = log)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.2488055 0.2309908 -5.406 6.43e-08 ***
## PA_invasive_presencewithout_invasive 0.4433465 0.1141518 3.884 0.000103 ***
## inside_PAoutside 0.0009267 0.0877093 0.011 0.991570
## NDVI 3.1490706 0.3680630 8.556 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(3.1445) family taken to be 1)
##
## Null deviance: 509.7 on 385 degrees of freedom
## Residual deviance: 385.8 on 382 degrees of freedom
## (4 observations deleted due to missingness)
## AIC: 1602.9
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 3.145
## Std. Err.: 0.465
##
## 2 x log-likelihood: -1592.857simms_ri <- simulateResiduals(fittedModel = glm_richness)
plotQQunif(simms_ri)
vif(glm_richness)## PA_invasive_presence inside_PA NDVI
## 1.123276 1.015078 1.124307tidy(glm_richness, conf.int = TRUE)## # A tibble: 4 × 7
## term estimate std.error statistic p.value conf.low conf.high
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) -1.25e+0 0.231 -5.41 6.43e- 8 -1.71 -0.797
## 2 PA_invasive_presence… 4.43e-1 0.114 3.88 1.03e- 4 0.231 0.658
## 3 inside_PAoutside 9.27e-4 0.0877 0.0106 9.92e- 1 -0.171 0.173
## 4 NDVI 3.15e+0 0.368 8.56 1.17e-17 2.45 3.85Plotting from Native Richness Predictions
gg_ndvi_rich <- ggeffect(glm_richness, terms = c("NDVI [all]"))
#Create NDVI plot with inside_PA using ggeffects
ggplot(gg_ndvi_rich, aes(x = x, y = predicted)) +
geom_line() +
geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = group), alpha = 0.2, colour = NA) +
theme_minimal() +
labs(x = "NDVI",y = "Predicted Native Species Richness",fill = "Confidence Intervals")
gg_invasive_rich <- ggeffect(glm_richness, terms = c("PA_invasive_presence [all]"))
# Create Invasive Abundance plot with inside_PA using ggeffects
ggplot(gg_invasive_rich, aes(x = x, y = predicted)) +
geom_bar(stat = "identity", position = "dodge", aes(fill = x)) +
geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2) +
theme_minimal() +
labs(x = "Invasive Rodent Presence", y = "Predicted Native Species Richness", fill = "Invasive Rodent Presence") +
scale_fill_manual(values = c("without_invasive" = "lightblue", "with_invasive" = "salmon")) 
# Print the predicted values along with their confidence intervals that are the error bars
print(gg_invasive_rich)## # Predicted counts of native_species_richness
##
## PA_invasive_presence | Predicted | 95% CI
## ---------------------------------------------
## with_invasive | 1.94 | 1.59, 2.35
## without_invasive | 3.02 | 2.73, 3.33Native diversity
GLM_sha<- lm(shannon_diversity ~ invasive_abundance + inside_PA + NDVI, data = trap_summ)
summary(GLM_sha)##
## Call:
## lm(formula = shannon_diversity ~ invasive_abundance + inside_PA +
## NDVI, data = trap_summ)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.05213 -0.44064 -0.06351 0.44794 1.48355
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.434416 0.163224 -2.661 0.00815 **
## invasive_abundance -0.009087 0.004541 -2.001 0.04619 *
## inside_PAoutside -0.051409 0.064150 -0.801 0.42347
## NDVI 1.801934 0.247892 7.269 2.56e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5712 on 336 degrees of freedom
## (50 observations deleted due to missingness)
## Multiple R-squared: 0.1684, Adjusted R-squared: 0.161
## F-statistic: 22.68 on 3 and 336 DF, p-value: 2.145e-13simms_sha <- simulateResiduals(fittedModel = GLM_sha)
plotQQunif(simms_sha)
vif(GLM_sha)## invasive_abundance inside_PA NDVI
## 1.061968 1.029379 1.060038tidy(GLM_sha, conf.int = TRUE)## # A tibble: 4 × 7
## term estimate std.error statistic p.value conf.low conf.high
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) -0.434 0.163 -2.66 8.15e- 3 -0.755 -0.113
## 2 invasive_abundance -0.00909 0.00454 -2.00 4.62e- 2 -0.0180 -0.000154
## 3 inside_PAoutside -0.0514 0.0641 -0.801 4.23e- 1 -0.178 0.0748
## 4 NDVI 1.80 0.248 7.27 2.56e-12 1.31 2.29Plotting from Native Rodent Diversity Predictions
# Generate predictions for NDVI
gg_ndvi_sha <- ggpredict(GLM_sha, terms = c("NDVI [all]"))
# Generate predictions for Invasive Abundance
gg_invasive_sha<- ggpredict(GLM_sha, terms = c("invasive_abundance [all]"))
# Plot for NDVI
sha_ndvi_plot <- ggplot(gg_ndvi_sha, aes(x = x, y = predicted)) +
geom_line() +
geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = group), alpha = 0.2, colour = NA) +
theme_minimal() +
labs(x = "NDVI", y = "Predicted Native Rodent Shannon Index", fill = "Confidence Intervals")
# Plot for Invasive Abundance
sha_invasive_plot <- ggplot(gg_invasive_sha, aes(x = x, y = predicted)) +
geom_line() +
geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = group), alpha = 0.2, colour = NA) +
theme_minimal() +
labs(x = "Invasive Abundance", y = "Predicted Native Rodent Shannon Index", fill = "Confidence Intervals")
# Combine plots
grid.arrange(sha_ndvi_plot, sha_invasive_plot, nrow = 2)
Environmental Factors Influencing Invasive Rodents
Invasive Rodent Abundance
Objective: Identify environmental factors associated with the relative abundance of invasive rodents, including NDVI and proximity to protected areas.
NDVI_mod <- glm.nb( invasive_abundance ~ NDVI + inside_PA, data = trap_summ)
summary(NDVI_mod)##
## Call:
## glm.nb(formula = invasive_abundance ~ NDVI + inside_PA, data = trap_summ,
## init.theta = 0.1286532515, link = log)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.9988 0.7450 8.052 8.16e-16 ***
## NDVI -8.4766 1.1957 -7.089 1.35e-12 ***
## inside_PAoutside -1.4339 0.3234 -4.434 9.26e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.1287) family taken to be 1)
##
## Null deviance: 284.25 on 385 degrees of freedom
## Residual deviance: 221.33 on 383 degrees of freedom
## (4 observations deleted due to missingness)
## AIC: 1035.2
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.1287
## Std. Err.: 0.0165
##
## 2 x log-likelihood: -1027.2140vif(NDVI_mod)## NDVI inside_PA
## 1.011206 1.011206tidy(NDVI_mod, conf.int = TRUE)## # A tibble: 3 × 7
## term estimate std.error statistic p.value conf.low conf.high
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 6.00 0.745 8.05 8.16e-16 4.57 7.51
## 2 NDVI -8.48 1.20 -7.09 1.35e-12 -10.8 -6.12
## 3 inside_PAoutside -1.43 0.323 -4.43 9.26e- 6 -2.06 -0.790Plotting from Invasive Abundance Predictions
gg_ndvi_inv <- ggeffect(NDVI_mod, terms = c("NDVI [all]", "inside_PA"))
#Create NDVI plot with inside_PA using ggeffects
ggplot(gg_ndvi_inv, aes(x = x, y = predicted, colour = group)) +
geom_line() +
geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = group), alpha = 0.2, colour = NA) +
theme_minimal() +
labs(x = "NDVI",y = "Predicted Invasive Abundance",colour = "Inside/Outside PA",fill = "Confidence Intervals") +
scale_y_continuous(expand = expansion(mult = c(0.1, 0.1)))
Tidy model outputs and create tables
# Tidy the model summaries for each of the models
tidy_glm_richness <- tidy(glm_richness)
tidy_NDVI_mod <- tidy(NDVI_mod)
tidy_GLM_sha <- tidy(GLM_sha)
tidy_presence_model <- tidy(presence_model)
tidy_abundance_model <- tidy(abundance_model)
# Create and save the table
tidy_presence_model %>%
kable("html", caption = "Native Presence Model") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = F) %>%
column_spec(1, bold = TRUE, color = "white", background = "lightblue") %>%
save_kable(file = "presence_model.html")
tidy_abundance_model %>%
kable("html", caption = "Native Abundance Model") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = F) %>%
column_spec(1, bold = TRUE, color = "white", background = "lightblue") %>%
save_kable(file = "abun_model.html")
tidy_glm_richness %>%
kable("html", caption = "Native Species Richness Model") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = F) %>%
column_spec(1, bold = TRUE, color = "white", background = "lightblue") %>%
save_kable(file = "rich_model.html")
tidy_GLM_sha %>%
kable("html", caption = "Native Shannon Index Model") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = F) %>%
column_spec(1, bold = TRUE, color = "white", background = "lightblue") %>%
save_kable(file = "shannon_model.html")
tidy_NDVI_mod %>%
kable("html", caption = "NDVI Model") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = F) %>%
column_spec(1, bold = TRUE, color = "white", background = "lightblue") %>%
save_kable(file = "NDVI_model.html")