This code calculates the geographic coverage (number of records per unit area) of the SCAR Antarctic Terrestrial Biodiversity Database. It also assesses the spatial autocorrelation of this coverage by calculating Moran’s I.

Initial setup

Load packages and set working directory.

library(tidyverse)

packages <- c("sf", "terra", "here", "viridis", "spdep", "ggpubr")

walk(packages, require, character.only = T)

here::here()

Load data

Read in the cleaned biodiversity database (cleaned in file “Terrestrial_Antarctic_Biodiversity_Database_Cleaning.Rmd”), the cleaned ice-free area vector (cleaned in file “ADD_ice_free_area_cleaning.Rmd”) and the coastline (Antarctic Digital Database coastline v. 7.5).

# Load biodiversity data as dataframe
Ant_biodf <- read.csv(here("Data/Ant_Terr_Bio_Data_Uncertainty_Terms_Removed_August_2023.csv"))

# Convert biodiversity dataframe to sf polygon
Ant_bio <- st_as_sf(Ant_biodf, 
                    coords = c("decimalLongitude", "decimalLatitude"), 
                    crs = 4326) # set as WGS 1984

# Project to WGS_1984_Stereographic_South_Pole
Ant_bio <- st_transform(Ant_bio, 3031) 

# Convert to a SpatVector for extract function below
Ant_bio_SPVE <- vect(Ant_bio)

# Load cleaned ice-free area vector 
ice_free <- st_read(here("Data/AAD_ice_free_vector_TIDIED.shp"), crs = 3031)

# Also convert to SpatVector for extract function
ice_free_SPVE <- vect(ice_free)

# Load the Antarctic coastline for plotting
coast <- st_read(here("Data/add_coastline_high_res_polygon_v7_5.shp"), crs = 3031)

# Convert to SpatVector for extract function
coastline_SPVE <- vect(coast)

# List target grid resolutions (1km, 20km, 100km, 200km)
res <- list(20000, 100000, 200000) 
# res <- list(1000, 20000, 100000, 200000) # 1km is slow to run and was run on QUT HPC. The following code is just for the 20km, 100km, and 200km resolutions.

# Index list
i <- list(1,2,3)
# i <- list(1,2,3,4)

Calculate geographic coverage

This function extracts the number of records found in each grid cell across resolutions. Adapted from code by Wauchope et al. (2019).

Coverage_func <- function(res, i) {

  # Count records per grid cell  ---------------------------------------------
    
# Create a 1km raster with all values as 1
grid <- rast(ext = ext(coastline_SPVE),
                   resolution = res, 
                   vals = 1) 

crs(grid) <- "epsg:3031"

# Add a one cell buffer to make sure all records are covered
grid <- terra::extend(grid, 1, fill = 1)

# Make raster and give all raster values a unique ID 
gridUNIQUE <- rast(ext = ext(grid), 
                   resolution = res, 
                   vals = 1:ncell(grid)) 

crs(gridUNIQUE) <- "epsg:3031"

# Find the ice-free grid cell ID each bio point lies in with a unique value
# ID is the ID of the ant bio observation, lyr.1 is the unique ID of the cell
extractUNIQUE <- terra::extract(gridUNIQUE, Ant_bio_SPVE)

# Count the number of points per unique grid cell
count <- count(extractUNIQUE, lyr.1)

# Make gridUNIQUE raster a data frame
gridUNIQUEdf <- as.data.frame(gridUNIQUE, xy = TRUE)

# Full join retains all values, join the grid dataframe and  the count values by the cell unique ID
update_rast <- full_join(gridUNIQUEdf, count, by = "lyr.1") %>% 
  .[, -3] %>% 
  drop_na(x) %>% # Remove if there's a row with no xy (where the data didn't overlap the coastline)
  replace(is.na(.), 0) %>% # Replace NAs (non-bio locs) with 0
  rast(.)  

crs(update_rast) <- "epsg:3031"

# Mask out non ice-free (ice-covered) areas ----------------------------------

# Set the grid cells that are NOT ice-free areas to 0
ice_free_mask <- terra::mask(grid, ice_free_SPVE, 
                      updatevalue = 0) 

# Mask the final raster of record counts by ice_free polygon, non ice-free areas become NA
update_rast_masked <- mask(update_rast, ice_free_mask,  
                           maskvalues = 0, updatevalue = NA)

coverage_rast_stack <<- c(coverage_rast_stack, update_rast_masked)

coverage_raster_list[[i]] <<- update_rast_masked

# Summarize total area covered vs. total area available -------------------

# Freq function counts the number of cells with x number of records: value is the number of records, count is the number of cells with that many records
# This no longer calculates including the ice-free cells as they are NA
freq <- as.data.frame(freq(update_rast_masked)) %>% 
  mutate(Perc_area_km2 = (count/sum(count))*100) %>% 
  mutate(Resolution = as.factor(res))

freq_all <<- dplyr::bind_rows(freq_all, freq)

}

# Run the function  ---------------------------------------------

# Set up empty objects to save results
coverage_raster_list <- vector(mode = "list", length = 4)
freq_all <- NULL
coverage_rast_stack <- c()

# Run function
map2(res, i , Coverage_func)

Optional - saving log of geographic coverage rasters for further plotting and plotting rasters in R.

# rast1 <- log(coverage_rast_stack[[1]] + 1)

rast2 <- log(coverage_rast_stack[[1]] + 1)

rast3 <- log(coverage_rast_stack[[2]] + 1)

rast4 <- log(coverage_rast_stack[[3]] + 1)
# 
# writeRaster(rast1, here("Outputs/Geographic_coverage_LOG_1km_rast.asc", overwrite = T))
# 
# writeRaster(rast2, here("Outputs/Geographic_coverage_LOG_20km_rast.asc", overwrite = T))
# 
# writeRaster(rast3, here("Outputs/Geographic_coverage_LOG_100km_rast.asc", overwrite = T))
# 
# writeRaster(rast4, here("Outputs/Geographic_coverage_LOG_200km_rast.asc", overwrite = T))


# Save each resolution as a dataframe
# geog_1kmdf <- as.data.frame(coverage_rast_stack[[1]], xy = TRUE)

geog_20kmdf <- as.data.frame(coverage_rast_stack[[1]], xy = TRUE)

geog_100kmdf <- as.data.frame(coverage_rast_stack[[2]], xy = TRUE)

geog_200kmdf <- as.data.frame(coverage_rast_stack[[3]], xy = TRUE)

# Set dataframe zeros to NA
# geog_1kmdf$n_nas <- ifelse(geog_1kmdf$n == 0, NA, geog_1kmdf$n)

geog_20kmdf$n_nas <- ifelse(geog_20kmdf$n == 0, NA, geog_20kmdf$n)

geog_100kmdf$n_nas <- ifelse(geog_100kmdf$n == 0, NA, geog_100kmdf$n)

geog_200kmdf$n_nas <- ifelse(geog_200kmdf$n == 0, NA, geog_200kmdf$n)

# Create plot
# plot1km <- ggplot() +
#   geom_sf(data = coast, fill = "grey94", alpha = 0.5, size = 0.00001)  +
#   geom_raster(data = geog_1kmdf, aes(x, y, fill = n_nas), alpha = 0.87) +
#   scale_fill_viridis(na.value = "black") +
#   labs(x = "Longitude", y = "Latitude", fill = str_wrap("Number of records", 17)) +
#   geom_area(inherit.aes = F, mapping = aes(y = 1, x = 1, color = "Zero records"), data = data.frame()) +
#   scale_colour_manual(name = element_blank(), values="#FFFFCC") +
#   guides(fill = guide_colourbar(ticks = F, order = 1), colour = guide_legend(order = 2)) +
#   theme_void()
# 
# plot1km <- plot1km + annotate("text", x = Inf, y = Inf, label = "1 km resolution", vjust = 1, hjust = 1)

plot20km <- ggplot() +
  geom_sf(data = coast, fill = "grey94", alpha = 0.5, size = 0.00001)  +
  geom_raster(data = geog_20kmdf, aes(x, y, fill = n_nas), alpha = 0.87) +
  scale_fill_viridis(na.value = "black") +
  labs(x = "Longitude", y = "Latitude", fill = str_wrap("Number of records", 17)) +
  geom_area(inherit.aes = F, mapping = aes(y = 1, x = 1, color = "Zero records"), data = data.frame()) +
  scale_colour_manual(name = element_blank(), values="#FFFFCC") +
  guides(fill = guide_colourbar(ticks = F, order = 1), colour = guide_legend(order = 2)) +
  theme_void()

plot20km <- plot20km + annotate("text", x = Inf, y = Inf, label = "20 km resolution", vjust = 1, hjust = 1)

plot100km <- ggplot() +
  geom_sf(data = coast, fill = "grey94", alpha = 0.5, size = 0.00001)  +
  geom_raster(data = geog_100kmdf, aes(x, y, fill = n_nas), alpha = 0.87) +
  scale_fill_viridis(na.value = "black") +
  labs(x = "Longitude", y = "Latitude", fill = str_wrap("Number of records", 17)) +
  geom_area(inherit.aes = F, mapping = aes(y = 1, x = 1, color = "Zero records"), data = data.frame()) +
  scale_colour_manual(name = element_blank(), values="#FFFFCC") +
  guides(fill = guide_colourbar(ticks = F, order = 1), colour = guide_legend(order = 2)) +
  theme_void()

plot100km <- plot100km + annotate("text", x = Inf, y = Inf, label = "100 km resolution", vjust = 1, hjust = 1)

plot200km <- ggplot() +
  geom_sf(data = coast, fill = "grey94", alpha = 0.5, size = 0.00001)  +
  geom_raster(data = geog_200kmdf, aes(x, y, fill = n_nas), alpha = 0.87) +
  scale_fill_viridis(na.value = "black") +
  labs(x = "Longitude", y = "Latitude", fill = str_wrap("Number of records", 17)) +
  geom_area(inherit.aes = F, mapping = aes(y = 1, x = 1, color = "Zero records"), data = data.frame()) +
  scale_colour_manual(name = element_blank(), values="#FFFFCC") +
  guides(fill = guide_colourbar(ticks = F, order = 1), colour = guide_legend(order = 2)) +
  theme_void()

plot200km <- plot200km + annotate("text", x = Inf, y = Inf, label = "200 km resolution", vjust = 1, hjust = 1)

# Arrange all plots into a figure
# geographic_coverage_plot <- ggarrange(plot1km + rremove("ylab") + rremove("xlab"), plot20km + rremove("ylab") + rremove("xlab"), plot100km + rremove("ylab") + rremove("xlab"), plot200km + rremove("ylab") + rremove("xlab"), nrow = 2, ncol = 2, align = "v", common.legend = TRUE, legend = "right") 

geographic_coverage_plot <- ggarrange(plot20km + rremove("ylab") + rremove("xlab"), plot100km + rremove("ylab") + rremove("xlab"), plot200km + rremove("ylab") + rremove("xlab"), nrow = 2, ncol = 2, align = "v", common.legend = TRUE, legend = "right") 

print(geographic_coverage_plot)

Summarise what percentage of cells have certain numbers of records in a table (used to estimate what percentage of cells have zero records across resolutions).

# Summarise counts by intervals
freq_table <- freq_all %>% 
  mutate(ints = cut(.$value, breaks = c(-1, 0, 1, 5, 10, 50, 100, 200, 500, 1000, 2000, 6000))) %>% 
  group_by(ints, Resolution) %>% 
  summarise(tot = sum(Perc_area_km2))
  
print(freq_table)
## # A tibble: 32 × 3
## # Groups:   ints [11]
##    ints   Resolution   tot
##    <fct>  <fct>      <dbl>
##  1 (-1,0] 20000      76.0 
##  2 (-1,0] 1e+05      45.3 
##  3 (-1,0] 2e+05      30.7 
##  4 (0,1]  20000       4.37
##  5 (0,1]  1e+05       7.73
##  6 (0,1]  2e+05       6.75
##  7 (1,5]  20000       6.77
##  8 (1,5]  1e+05      10.7 
##  9 (1,5]  2e+05      14.1 
## 10 (5,10] 20000       3.18
## # ℹ 22 more rows

Spatial autocorrelation of coverage

Using Moran’s I from the package spdep, calculate the spatial autocorrelation of the grid cell coverage. First we set up the neighbourhood distance to define cells that are considered ‘adjacent’. We chose neighbourhood distances for our 1km, 20km, 100km, and 200km resolutions to be 1.5km, 30km, 150km, and 300km. Neighbours are therefore adjacent by the Queen’s adjacency rule.

# Set up the dataframe for saving
MoransI <- data.frame(Resolution = numeric(), MoransI = numeric())

MoransI_func <- function(res, i) {
  
  # Calculating Moran's I  -------------------------------------------------
  
  # Calculate the neighbourhood distance as 1.5 x resolution
  neighbour_dist <- res*1.5
  
  update_rast_xy <- coverage_raster_list[[i]] %>% 
    as.data.frame(xy=TRUE) %>% 
    .[, -3] #remove n
  
  print("update_rast_xy created")
  
  # dnearneigh is in metres, finds to index of all values within the specified distance
  nb1 <- dnearneigh(as.matrix(update_rast_xy), 0, neighbour_dist)
  
  # For viewing which neighbours are selected
  # plot(nb1, coords = as.matrix(update_rast_xy))
  
  print("...Nearest neighbours found...")
  
  # Print list of distances
  # dists <- nbdists(nb1, update_rast_xy)
  
  # Generate a weights object
  weights <- nb2listw(nb1, 
                      glist = NULL,
                      style = "B",
                      zero.policy = TRUE) 
  
  print("...Weights calculated...")
  
  # Then extract just the values per cell
  var <- coverage_raster_list[[i]] %>% 
    as.data.frame() %>% 
    .[,1] 
  
  # Normalize the data so it is not affected by high values
  normalize <- function(x, na.rm = TRUE) {
    return((x- min(x)) /(max(x)-min(x)))
  }
  
  norm_var <- normalize(var)
  
  # Moran's I value for normalized data
  I <- moran(norm_var, 
             weights, 
             zero.policy = T, 
             length(nb1), 
             Szero(weights)) %>% 
    .[[1]]
  
  print("...Moran's I calculated...")
    # Add output to dataframe
  MoransI <<- MoransI %>% 
    add_row(Resolution = res, MoransI = I)
  
  # Hypothesis test  ----------------------------------------------------
  
  # Run a hypothesis test via permutation
  set.seed(1)
  moran_bootstrap <- moran.mc(var,
                              weights,
                              nsim = 10000,
                              zero.policy = TRUE)
  
  print("...Bootstrapping finished...")
  
  # Print summary outputs  ----------------------------------------------
  
  # Print the output of the hypothesis test
  print(moran_bootstrap)
  
  # Plot the density of expected (under randomisation) vs. observed values of Moran's I 
  plot(moran_bootstrap, main="", las=1)
  
  # Plot the relationship between a cell's value (no. records) and its nearest neighbour values
  moran.plot(var, weights, zero.policy = T)
  
  # Calculate the number of bins for a correlogram according to Sturge's rule (refer to Borcard et al. 2018)
  no_recs <- coverage_raster_list[[i]] %>% 
    as.data.frame() %>% 
    pull(n) %>% 
    sum 
  
  no_dists <- round(1 + 3.322*log10(no_recs))
  
  # Generate a correlogram to examine the spatial scale of correlation
  subs.correlog <-
    sp.correlogram(nb1, 
                   norm_var, 
                   order = no_dists, 
                   method = "I", 
                   zero.policy = TRUE)
  
  # Print the results with a "Holm" correction for multiple testing
  print(subs.correlog, p.adj.method = "bonferroni")
  
  plot(subs.correlog)
  
}

map2(res, i, MoransI_func)
## [1] "update_rast_xy created"
## [1] "...Nearest neighbours found..."
## [1] "...Weights calculated..."
## [1] "...Moran's I calculated..."
## [1] "...Bootstrapping finished..."
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  var 
## weights: weights  
## number of simulations + 1: 10001 
## 
## statistic = 0.24986, observed rank = 10001, p-value = 9.999e-05
## alternative hypothesis: greater

## Spatial correlogram for norm_var 
## method: Moran's I
##              estimate expectation    variance standard deviate Pr(I) two sided
## 1 (3169)   2.2734e-01 -3.1566e-04  1.2637e-04          20.2516       < 2.2e-16
## 2 (3075)   1.4586e-01 -3.2531e-04  1.0037e-04          14.5917       < 2.2e-16
## 3 (2953)   1.1066e-01 -3.3875e-04  8.9756e-05          11.7160       < 2.2e-16
## 4 (2876)   8.6881e-02 -3.4783e-04  8.6481e-05           9.3799       < 2.2e-16
## 5 (2820)   7.5895e-02 -3.5474e-04  8.6300e-05           8.2079       3.602e-15
## 6 (2779)   9.2170e-02 -3.5997e-04  8.7764e-05           9.8770       < 2.2e-16
## 7 (2728)   9.3766e-02 -3.6670e-04  9.1243e-05           9.8546       < 2.2e-16
## 8 (2665)   5.1044e-02 -3.7538e-04  9.3694e-05           5.3121       1.734e-06
## 9 (2588)   1.7126e-02 -3.8655e-04  9.3989e-05           1.8064        1.000000
## 10 (2501)  2.2604e-03 -4.0000e-04  9.3113e-05           0.2757        1.000000
## 11 (2412)  5.1132e-03 -4.1477e-04  8.7350e-05           0.5915        1.000000
## 12 (2333)  5.3460e-03 -4.2882e-04  8.3197e-05           0.6331        1.000000
## 13 (2276)  1.4396e-02 -4.3956e-04  8.0693e-05           1.6515        1.000000
## 14 (2228)  1.8697e-02 -4.4904e-04  7.9505e-05           2.1472        0.508431
## 15 (2187)  1.7424e-02 -4.5746e-04  7.8799e-05           2.0144        0.703531
## 16 (2155)  3.2830e-02 -4.6425e-04  8.1233e-05           3.6941        0.003531
##              
## 1 (3169)  ***
## 2 (3075)  ***
## 3 (2953)  ***
## 4 (2876)  ***
## 5 (2820)  ***
## 6 (2779)  ***
## 7 (2728)  ***
## 8 (2665)  ***
## 9 (2588)     
## 10 (2501)    
## 11 (2412)    
## 12 (2333)    
## 13 (2276)    
## 14 (2228)    
## 15 (2187)    
## 16 (2155) ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## [1] "update_rast_xy created"
## [1] "...Nearest neighbours found..."
## [1] "...Weights calculated..."
## [1] "...Moran's I calculated..."
## [1] "...Bootstrapping finished..."
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  var 
## weights: weights  
## number of simulations + 1: 10001 
## 
## statistic = 0.24172, observed rank = 10000, p-value = 9.999e-05
## alternative hypothesis: greater

## Spatial correlogram for norm_var 
## method: Moran's I
##            estimate expectation   variance standard deviate Pr(I) two sided    
## 1 (370)   0.2166767  -0.0027100  0.0011013           6.6109       6.111e-10 ***
## 2 (362)   0.1655489  -0.0027701  0.0010100           5.2964       1.890e-06 ***
## 3 (354)   0.2007731  -0.0028329  0.0010493           6.2854       5.232e-09 ***
## 4 (354)   0.1576367  -0.0028329  0.0011457           4.7409       3.405e-05 ***
## 5 (349)   0.0437933  -0.0028736  0.0012652           1.3120               1    
## 6 (340)   0.0482899  -0.0029499  0.0014203           1.3596               1    
## 7 (327)   0.0206006  -0.0030675  0.0014970           0.6117               1    
## 8 (312)  -0.0301507  -0.0032154  0.0015763          -0.6784               1    
## 9 (295)  -0.0502347  -0.0034014  0.0017354          -1.1242               1    
## 10 (275) -0.0585699  -0.0036496  0.0018355          -1.2819               1    
## 11 (247) -0.0676104  -0.0040650  0.0020692          -1.3969               1    
## 12 (216) -0.0489380  -0.0046512  0.0024460          -0.8955               1    
## 13 (189) -0.0343710  -0.0053191  0.0028518          -0.5440               1    
## 14 (167) -0.0180347  -0.0060241  0.0032741          -0.2099               1    
## 15 (144)  0.0060722  -0.0069930  0.0036528           0.2162               1    
## 16 (119) -0.0027471  -0.0084746  0.0038717           0.0920               1    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## [1] "update_rast_xy created"
## [1] "...Nearest neighbours found..."
## [1] "...Weights calculated..."
## [1] "...Moran's I calculated..."
## [1] "...Bootstrapping finished..."
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  var 
## weights: weights  
## number of simulations + 1: 10001 
## 
## statistic = 0.39437, observed rank = 10001, p-value = 9.999e-05
## alternative hypothesis: greater

## Spatial correlogram for norm_var 
## method: Moran's I
##            estimate expectation   variance standard deviate Pr(I) two sided    
## 1 (162)   0.3745894  -0.0062112  0.0025453           7.5480       7.075e-13 ***
## 2 (160)   0.3846369  -0.0062893  0.0023745           8.0225       1.658e-14 ***
## 3 (160)   0.0887519  -0.0062893  0.0025986           1.8644          0.9962    
## 4 (157)   0.0187631  -0.0064103  0.0025434           0.4992          1.0000    
## 5 (155)  -0.0395626  -0.0064935  0.0024980          -0.6616          1.0000    
## 6 (153)  -0.0385488  -0.0065789  0.0025176          -0.6372          1.0000    
## 7 (153)  -0.0323806  -0.0065789  0.0024601          -0.5202          1.0000    
## 8 (153)  -0.0303159  -0.0065789  0.0026071          -0.4649          1.0000    
## 9 (153)  -0.0488184  -0.0065789  0.0026889          -0.8146          1.0000    
## 10 (151) -0.0591356  -0.0066667  0.0027098          -1.0079          1.0000    
## 11 (148) -0.0704868  -0.0068027  0.0026795          -1.2303          1.0000    
## 12 (145) -0.0833541  -0.0069444  0.0027980          -1.4445          1.0000    
## 13 (141) -0.0745712  -0.0071429  0.0030864          -1.2137          1.0000    
## 14 (136) -0.0507531  -0.0074074  0.0033466          -0.7493          1.0000    
## 15 (130)  0.0064579  -0.0077519  0.0036432           0.2354          1.0000    
## 16 (120)  0.0913509  -0.0084034  0.0040731           1.5630          1.0000    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## [[1]]
## NULL
## 
## [[2]]
## NULL
## 
## [[3]]
## NULL
# Print the dataframe of estimates
print(MoransI)
##   Resolution   MoransI
## 1      2e+04 0.2525409
## 2      1e+05 0.2449882
## 3      2e+05 0.3968023

References

Borcard, D., Gillet, F. and Legendre, P., 2018. Numerical ecology with R (Vol. 2, p. 688). New York: Springer.

Wauchope, H.S., Shaw, J.D., Terauds, A., 2019. A snapshot of biodiversity protection in Antarctica. Nature Communications 10, 946. https://doi.org/10.1038/s41467-019-08915-6

Session information

sessionInfo()
## R version 4.2.3 (2023-03-15 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19045)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=English_Australia.utf8  LC_CTYPE=English_Australia.utf8   
## [3] LC_MONETARY=English_Australia.utf8 LC_NUMERIC=C                      
## [5] LC_TIME=English_Australia.utf8    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] ggpubr_0.6.0      spdep_1.2-8       spData_2.2.2      viridis_0.6.2    
##  [5] viridisLite_0.4.1 here_1.0.1        terra_1.7-18      sf_1.0-12        
##  [9] lubridate_1.9.2   forcats_1.0.0     stringr_1.5.0     dplyr_1.1.1      
## [13] purrr_1.0.1       readr_2.1.4       tidyr_1.3.0       tibble_3.2.1     
## [17] ggplot2_3.4.2     tidyverse_2.0.0  
## 
## loaded via a namespace (and not attached):
##  [1] sass_0.4.5         jsonlite_1.8.4     carData_3.0-5      bslib_0.4.2       
##  [5] highr_0.10         sp_1.6-0           yaml_2.3.7         pillar_1.9.0      
##  [9] backports_1.4.1    lattice_0.20-45    glue_1.6.2         digest_0.6.31     
## [13] ggsignif_0.6.4     colorspace_2.1-0   cowplot_1.1.1      htmltools_0.5.5   
## [17] pkgconfig_2.0.3    broom_1.0.4        s2_1.1.2           scales_1.2.1      
## [21] tzdb_0.3.0         timechange_0.2.0   proxy_0.4-27       generics_0.1.3    
## [25] farver_2.1.1       car_3.1-2          cachem_1.0.7       withr_2.5.0       
## [29] cli_3.6.1          magrittr_2.0.3     deldir_1.0-6       evaluate_0.20     
## [33] fansi_1.0.4        rstatix_0.7.2      class_7.3-21       tools_4.2.3       
## [37] hms_1.1.3          lifecycle_1.0.3    munsell_0.5.0      compiler_4.2.3    
## [41] jquerylib_0.1.4    e1071_1.7-13       rlang_1.1.0        classInt_0.4-9    
## [45] units_0.8-1        grid_4.2.3         rstudioapi_0.14    labeling_0.4.2    
## [49] rmarkdown_2.21     boot_1.3-28.1      wk_0.7.2           gtable_0.3.3      
## [53] codetools_0.2-19   abind_1.4-5        DBI_1.1.3          R6_2.5.1          
## [57] gridExtra_2.3      knitr_1.42         fastmap_1.1.1      utf8_1.2.3        
## [61] rprojroot_2.0.3    KernSmooth_2.23-20 stringi_1.7.12     Rcpp_1.0.10       
## [65] vctrs_0.6.1        tidyselect_1.2.0   xfun_0.38