Hourly summaries of the observed and simulated trajecotries
Load packages
Import the covariate raster layers that were used in the step selection model fitting.
Code
# unscaled rasters
ndvi_stack <- rast("mapping/cropped rasters/ndvi_GEE_projected_watermask20230207.tif")
ndvi_2018_dry <- ndvi_stack[[8:10]]
ndvi_2019_dry <- ndvi_stack[[19:21]]
ndvi_dry <- terra::mean(c(ndvi_2018_dry, ndvi_2019_dry))
names(ndvi_dry) <- "ndvi_dry"
# plot(ndvi_dry)
canopy <- rast("mapping/cropped rasters/canopy_cover.tif")
herby <- rast("mapping/cropped rasters/veg_herby.tif")
slope <- rast("mapping/cropped rasters/slope_raster.tif")Load simulated trajectories and extract covariate information
In this script we are only calculating the summaries for the simulations of a single model. To compare the simulations of the other models, we would just change the folder name that we are reading in the csv files from.
We have the trajectories saved in certain folders, so we make a list of all the files in that folder, and then filter them based on the trajectories that we want to select (using the grep function to match the strings). Each file will have many trajectories for the starting locations of each of the buffalo (depending on how the simulation was set up). We then read in the csv files.
Code
# read in multiple csv files with similar filenames and bind them together
sim_data_full_list <-
list.files("Rev1_simulations/no_memory_fit/3p/",
pattern = "*.csv", full.names = T)
sim_data_all <- grep("50ch", sim_data_full_list[1:10], value = T) %>%
grep("3000", x = ., value = T) %>%
map_dfr(read_csv, .id = "file_num") %>%
mutate(file_num = str_remove(file_num, ".csv")) %>%
mutate(file_num = as.numeric(file_num))
# to check the initial locations of the simulated trajectories
sim_data_all %>% group_by(id) %>% slice(1)Prepare the simulated data
We need to do a bit of tidying, mostly just so we have a unique ID for each of the simulated trajectories, and then we convert to a track object using the amt package and extract covariate values.
Code
# create vector of unique ids for subsetting
sim_data_ids <- unique(sim_data_all$id)
# attach ACST time format
tz(sim_data_all$t_) <- "Australia/Darwin"
# subtract a year from the time column so it matches the buffalo data
# minus the warmup period (500 hours)
sim_data_all$t_ <- sim_data_all$t_ - years(1) - days(31)
# convert to track object for step lengths and turning angles etc
sim_data_all <- sim_data_all %>% mk_track(id = id, x_, y_, t_, all_cols = T, crs = 3112)
sim_data_all_nested <- sim_data_all %>% arrange(id) %>% nest(data = -"id")
# extract covariate values at the end of the all the simulated steps
sim_data_all_nested_steps <- sim_data_all_nested %>%
mutate(steps = map(data, function(x)
x %>% steps(keep_cols = "end") %>%
extract_covariates(ndvi_dry,
where = "end") %>%
extract_covariates(canopy,
where = "end") %>%
extract_covariates(herby,
where = "end") %>%
extract_covariates(slope,
where = "end")))
# sim_data_all_nested_steps$steps[[1]]
sim_data_all_steps <- sim_data_all_nested_steps %>%
amt::select(id, steps) %>%
amt::unnest(cols = steps)
# rearrange the columns so it is id num, sim num, file num then the rest
sim_data_all_steps <- sim_data_all_steps %>%
# dplyr::select(id, file_num, sim_num, traj, everything()) %>%
mutate(month = month(t1_))
sim_data_all_track <- sim_data_all_steps %>%
mk_track(id = id, x1_, y1_, t1_, order_by_ts = T, all_cols = T, crs = 3112) %>%
arrange(id)
sim_ids <- unique(sim_data_all_steps$id)
# unique(sim_data_all_steps$file_num)Split the track into bursts when it crosses to the other side of the landscape
As we had a wrapped landscape in the simulations, we need to split the trajectories into bursts when they go outside the extent and jump to the other side. We also remove any sections that have less than 3 points, and any steps that jump to the other side of the map (determined by steps greater than 10,000m, which is much larger than the largest observed step length).
Code
# Identify jumps by checking if the distance is greater than the threshold
jumps <- which(sim_data_all_steps$sl_ > 5000)
jumps <- c(1, jumps)
bursts <- rep(NA, nrow(sim_data_all_steps))
for(i in 1:(length(jumps)-1)) {
bursts[jumps[i]:jumps[i+1]] <- i
}
sim_data_all_steps <- sim_data_all_steps %>% mutate(burst = bursts)
# remove any sections that have less than 3 points
sim_data_all_steps <- sim_data_all_steps %>%
group_by(burst) %>%
filter(n() > 2) %>%
ungroup()
# remove steps that jump to the other side of the map
sim_data_all_steps <- sim_data_all_steps %>% filter(sl_ < 10000)
burst_ids <- unique(sim_data_all_steps$burst)
sim_data_all_track <- sim_data_all_steps %>%
mk_track(id = id, x1_, y1_, t1_, order_by_ts = T, all_cols = T, crs = 3112) %>%
arrange(id)Plot some example trajectories (by burst)
Load observed buffalo data for comparison
Code
Rows: 1165406 Columns: 22
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
dbl (18): id, burst_, x1_, x2_, y1_, y2_, sl_, ta_, dt_, hour_t2, step_id_,...
lgl (1): case_
dttm (3): t1_, t2_, t2_rounded
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.Code
buffalo_data_all <- buffalo_data_all %>%
mutate(t1_ = lubridate::with_tz(buffalo_data_all$t1_, tzone = "Australia/Darwin"),
t2_ = lubridate::with_tz(buffalo_data_all$t2_, tzone = "Australia/Darwin"))
buffalo_data_all <- buffalo_data_all %>%
mutate(id_num = as.numeric(factor(id)),
step_id = step_id_,
x1 = x1_, x2 = x2_,
y1 = y1_, y2 = y2_,
t1 = t1_,
t1_rounded = round_date(t1_, "hour"),
hour_t1 = hour(t1_rounded),
t2 = t2_,
t2_rounded = round_date(t2_, "hour"),
hour_t2 = hour(t2_rounded),
hour = ifelse(hour_t2 == 0, 24, hour_t2),
yday = yday(t1_),
year = year(t1_),
month = month(t1_),
sl = sl_,
log_sl = log(sl_),
ta = ta_,
cos_ta = cos(ta_),
canopy_01 = canopy_cover/100)
# ensure that we have just the individuals that were used for model fitting
buffalo_year_ids <- c(2005, 2014, 2018, 2021, 2022, 2024, 2039,
2154, 2158, 2223, 2327, 2387, 2393)
buffalo_data_all <- buffalo_data_all %>% filter(id %in% buffalo_year_ids)
buffalo_data_all <- buffalo_data_all %>% filter(t1 < "2019-07-25 09:32:42 ACST")
buffalo_ids <- unique(buffalo_data_all$id)
# check the initial location and time for each individual
buffalo_data_all %>% group_by(id) %>% slice(1)Code
# make a dataframe of only the presence locations
buffalo_data_all_pres <- buffalo_data_all %>% filter(y == 1)
# convert to track object for step lengths and turning angles etc
buffalo_data_all_pres <- buffalo_data_all_pres %>%
mk_track(id = id, x1_, y1_, t1_, order_by_ts = T, all_cols = T, crs = 3112) %>%
arrange(id) %>%
extract_covariates(ndvi_dry, where = "end")Subset the buffalo data to the period that simulations were generated for
Keep the simulated data that lines up with the observed data, for each individual buffalo
As we want the simulated trajectories to cover the same period as the observed data, we filter the simulated data between the minimum and maximum observed times for the buffalo
Hourly movement behaviour and selection of covariates
Here we bin the trajectories into the hours of the day, and calculate the mean, median (where appropriate) and sd values for the step lengths and four habitat covariates.
We save the results as a csv to compare between all of the models.
Observed buffalo data
Code
buffalo_hourly_habitat <-
buffalo_data %>% dplyr::group_by(hour, id) %>%
summarise(n = n(),
step_length_mean = mean(sl_),
step_length_median = median(sl_),
step_length_sd = sd(sl_),
ndvi_mean = mean(ndvi_dry),
ndvi_median = median(ndvi_dry),
ndvi_sd = sd(ndvi_dry),
herby_mean = mean(veg_herby),
herby_sd = sd(veg_herby),
canopy_mean = mean(canopy_cover/100),
canopy_sd = sd(canopy_cover/100),
slope_mean = mean(slope),
slope_median = median(slope),
slope_sd = sd(slope)
) %>% ungroup()`summarise()` has grouped output by 'hour'. You can override using the
`.groups` argument.Code
buffalo_hourly_habitat <- data.frame("data" = "buffalo", buffalo_hourly_habitat) %>%
mutate(id = as.factor(id))
# write the summaries to a csv
# write.csv(buffalo_hourly_habitat,
# paste0("outputs/buffalo_summaries_hourly_habitat.csv"))
buffalo_hourly_habitat_long <- buffalo_hourly_habitat %>%
pivot_longer(cols = !c(data, hour, id), names_to = "variable", values_to = "value")
sim_hourly_habitat <-
sim_data_btime %>% dplyr::group_by(hour, id) %>%
summarise(n = n(),
step_length_mean = mean(sl_),
step_length_median = median(sl_),
step_length_sd = sd(sl_),
ndvi_mean = mean(ndvi_dry),
ndvi_median = median(ndvi_dry),
ndvi_sd = sd(ndvi_dry),
herby_mean = mean(veg_herby),
herby_sd = sd(veg_herby),
canopy_mean = mean(canopy_cover/100),
canopy_sd = sd(canopy_cover/100),
slope_mean = mean(slope_raster),
slope_median = median(slope_raster),
slope_sd = sd(slope_raster)
) %>% ungroup()`summarise()` has grouped output by 'hour'. You can override using the
`.groups` argument.Code
# change the model name accordingly
sim_hourly_habitat <- data.frame("data" = "3p", sim_hourly_habitat) %>%
rename(id = id) %>% mutate(id = as.factor(id))
# save the summaries - change the model name accordingly
# write.csv(sim_hourly_habitat,
# paste0("outputs/sim_3p_summaries_hourly_habitat.csv"))
sim_hourly_habitat_long <- sim_hourly_habitat %>%
pivot_longer(cols = !c(data, hour, id), names_to = "variable", values_to = "value")
# combine the dataframe
hourly_habitat_long <- bind_rows(buffalo_hourly_habitat_long, sim_hourly_habitat_long)Comparing hourly summary statistics
Import the summary statistic data frames from the observed and simulated datasets that were outputted by the trajectory validation scripts, and combine into a single data frame for plotting
Code
# read in the observed data
summaries_hourly_buffalo <-
read_csv("outputs/buffalo_summaries_hourly_habitat.csv") %>%
mutate(id = as.factor(id))
# read in the simulated data
summaries_hourly_0p <-
read_csv("outputs/sim_0p_summaries_hourly_habitat.csv") %>%
mutate(id = as.factor(id))
summaries_hourly_1p <-
read_csv("outputs/sim_1p_summaries_hourly_habitat.csv") %>%
mutate(id = as.factor(id))
summaries_hourly_2p <-
read_csv("outputs/sim_2p_summaries_hourly_habitat.csv") %>%
mutate(id = as.factor(id))
summaries_hourly_3p <-
read_csv("outputs/sim_3p_summaries_hourly_habitat.csv") %>%
mutate(id = as.factor(id))
summaries_hourly_all <- bind_rows(summaries_hourly_buffalo,
summaries_hourly_0p,
summaries_hourly_1p,
summaries_hourly_2p,
summaries_hourly_3p)
summaries_hourly_all_long <- summaries_hourly_all %>%
dplyr::select(!...1) %>%
pivot_longer(cols = !c(data, hour, id),
names_to = "variable",
values_to = "value") %>%
mutate(Data = factor(str_to_title(data),
levels = c("Buffalo", "0p", "1p", "2p", "3p")))
head(summaries_hourly_all)Plot the hourly habitat selection between the observed and simulated data from different models
In these plots we are assessing how well the simulated trajectories capture the movement dynamics and habitat use of the observed data.
To express the stochasticity of the simulations, here we show the 25th to 50th quantiles and the 2.5th to 97.5th quantiles of the data. Remember that the ‘data’ are the means for each hour for each trajectory, so the quantiles are calculated across the means for each hour.
Here we create a ribbon from the 25th to 50th quantiles and from the 2.5th to 97.5th quantiles of the data
Code
hourly_summary_quantiles <- summaries_hourly_all_long %>%
dplyr::group_by(Data, hour, variable) %>%
summarise(n = n(),
mean = mean(value, na.rm = TRUE),
sd = sd(value, na.rm = TRUE),
q025 = quantile(value, probs = 0.025, na.rm = TRUE),
q25 = quantile(value, probs = 0.25, na.rm = TRUE),
q50 = quantile(value, probs = 0.5, na.rm = TRUE),
q75 = quantile(value, probs = 0.75, na.rm = TRUE),
q975 = quantile(value, probs = 0.975, na.rm = TRUE))`summarise()` has grouped output by 'Data', 'hour'. You can override using the
`.groups` argument.Here we show the hourly step length and selection of three of the habitat covariates. We use a dashed-line ribbon for the 95% interval and a solid-line ribbon for the 50% interval. We show the mean as a solid line.
Set plotting parameters
Code
# set plotting parameters here that will change in each plot
buff_path_alpha <- 0.1
ribbon_95_alpha <- 0.1
ribbon_50_alpha <- 0.15
path_95_alpha <- 1
# set path alpha
buff_path_alpha <- 0.25
# linewidth
buff_path_linewidth <- 0.5
# Create color mapping
unique_groups <- unique(summaries_hourly_all_long$Data)
colors <- viridis(length(unique_groups))
names(colors) <- unique_groups
colors["Buffalo"] <- "red"Code
hourly_path_sl_plot <- ggplot() +
geom_ribbon(data = hourly_summary_quantiles %>%
filter(Data == "Buffalo" & variable == "step_length_mean"),
aes(x = hour, ymin = q25, ymax = q75, fill = Data),
alpha = ribbon_50_alpha) +
geom_ribbon(data = hourly_summary_quantiles %>%
filter(!Data == "Buffalo" & variable == "step_length_mean"),
aes(x = hour, ymin = q25, ymax = q75, fill = Data),
alpha = ribbon_50_alpha) +
geom_path(data = summaries_hourly_all_long %>%
filter(Data == "Buffalo" & variable == "step_length_mean"),
aes(x = hour, y = value, colour = Data, group = interaction(id, Data)),
alpha = buff_path_alpha,
linewidth = buff_path_linewidth) +
geom_path(data = hourly_summary_quantiles %>%
filter(!Data == "Buffalo" & variable == "step_length_mean"),
aes(x = hour, y = q025, colour = Data),
linetype = "dashed",
alpha = path_95_alpha) +
geom_path(data = hourly_summary_quantiles %>%
filter(!Data == "Buffalo" & variable == "step_length_mean"),
aes(x = hour, y = q975, colour = Data),
linetype = "dashed",
alpha = path_95_alpha) +
geom_path(data = hourly_summary_quantiles %>%
filter(Data == "Buffalo" & variable == "step_length_mean"),
aes(x = hour, y = mean, colour = Data),
linewidth = 1) +
geom_path(data = hourly_summary_quantiles %>%
filter(!Data == "Buffalo" & variable == "step_length_mean"),
aes(x = hour, y = mean, colour = Data),
linewidth = 1) +
scale_fill_manual(values = colors) +
scale_colour_manual(values = colors) +
scale_x_discrete("Hour", breaks = seq(0,24,3)) +
scale_y_continuous("Mean value") +
ggtitle("Step length (m)") +
theme_classic()
hourly_path_sl_plot
Code
hourly_path_ndvi_plot <- ggplot() +
geom_ribbon(data = hourly_summary_quantiles %>%
filter(Data == "Buffalo" & variable == "ndvi_mean"),
aes(x = hour, ymin = q25, ymax = q75, fill = Data),
alpha = ribbon_50_alpha) +
geom_ribbon(data = hourly_summary_quantiles %>%
filter(!Data == "Buffalo" & variable == "ndvi_mean"),
aes(x = hour, ymin = q25, ymax = q75, fill = Data),
alpha = ribbon_50_alpha) +
geom_path(data = summaries_hourly_all_long %>%
filter(Data == "Buffalo" & variable == "ndvi_mean"),
aes(x = hour, y = value, colour = Data, group = interaction(id, Data)),
alpha = buff_path_alpha,
linewidth = buff_path_linewidth) +
geom_path(data = hourly_summary_quantiles %>%
filter(!Data == "Buffalo" & variable == "ndvi_mean"),
aes(x = hour, y = q025, colour = Data),
linetype = "dashed",
alpha = path_95_alpha) +
geom_path(data = hourly_summary_quantiles %>%
filter(!Data == "Buffalo" & variable == "ndvi_mean"),
aes(x = hour, y = q975, colour = Data),
linetype = "dashed",
alpha = path_95_alpha) +
geom_path(data = hourly_summary_quantiles %>%
filter(Data == "Buffalo" & variable == "ndvi_mean"),
aes(x = hour, y = mean, colour = Data),
linewidth = 1) +
geom_path(data = hourly_summary_quantiles %>%
filter(!Data == "Buffalo" & variable == "ndvi_mean"),
aes(x = hour, y = mean, colour = Data),
linewidth = 1) +
scale_fill_manual(values = colors) +
scale_colour_manual(values = colors) +
scale_x_discrete("Hour", breaks = seq(0,24,3)) +
scale_y_continuous("Mean value") +
ggtitle("NDVI") +
theme_classic()
hourly_path_ndvi_plot
Code
hourly_path_herby_plot <- ggplot() +
geom_ribbon(data = hourly_summary_quantiles %>%
filter(Data == "Buffalo" & variable == "herby_mean"),
aes(x = hour, ymin = q25, ymax = q75, fill = Data),
alpha = ribbon_50_alpha) +
geom_ribbon(data = hourly_summary_quantiles %>%
filter(!Data == "Buffalo" & variable == "herby_mean"),
aes(x = hour, ymin = q25, ymax = q75, fill = Data),
alpha = ribbon_50_alpha) +
geom_path(data = summaries_hourly_all_long %>%
filter(Data == "Buffalo" & variable == "herby_mean"),
aes(x = hour, y = value, colour = Data, group = interaction(id, Data)),
alpha = buff_path_alpha,
linewidth = buff_path_linewidth) +
geom_path(data = hourly_summary_quantiles %>%
filter(!Data == "Buffalo" & variable == "herby_mean"),
aes(x = hour, y = q025, colour = Data),
linetype = "dashed",
alpha = path_95_alpha) +
geom_path(data = hourly_summary_quantiles %>%
filter(!Data == "Buffalo" & variable == "herby_mean"),
aes(x = hour, y = q975, colour = Data),
linetype = "dashed",
alpha = path_95_alpha) +
geom_path(data = hourly_summary_quantiles %>%
filter(Data == "Buffalo" & variable == "herby_mean"),
aes(x = hour, y = mean, colour = Data),
linewidth = 1) +
geom_path(data = hourly_summary_quantiles %>%
filter(!Data == "Buffalo" & variable == "herby_mean"),
aes(x = hour, y = mean, colour = Data),
linewidth = 1) +
scale_fill_manual(values = colors) +
scale_colour_manual(values = colors) +
scale_x_discrete("Hour", breaks = seq(0,24,3)) +
scale_y_continuous("Mean value") +
ggtitle("Herbaceous vegetation") +
theme_classic()
hourly_path_herby_plot
Code
hourly_path_canopy_plot <- ggplot() +
geom_ribbon(data = hourly_summary_quantiles %>%
filter(Data == "Buffalo" & variable == "canopy_mean"),
aes(x = hour, ymin = q25, ymax = q75, fill = Data),
alpha = ribbon_50_alpha) +
geom_ribbon(data = hourly_summary_quantiles %>%
filter(!Data == "Buffalo" & variable == "canopy_mean"),
aes(x = hour, ymin = q25, ymax = q75, fill = Data),
alpha = ribbon_50_alpha) +
geom_path(data = summaries_hourly_all_long %>%
filter(Data == "Buffalo" & variable == "canopy_mean"),
aes(x = hour, y = value, colour = Data, group = interaction(id, Data)),
alpha = buff_path_alpha,
linewidth = buff_path_linewidth) +
geom_path(data = hourly_summary_quantiles %>%
filter(!Data == "Buffalo" & variable == "canopy_mean"),
aes(x = hour, y = q025, colour = Data),
linetype = "dashed",
alpha = path_95_alpha) +
geom_path(data = hourly_summary_quantiles %>%
filter(!Data == "Buffalo" & variable == "canopy_mean"),
aes(x = hour, y = q975, colour = Data),
linetype = "dashed",
alpha = path_95_alpha) +
geom_path(data = hourly_summary_quantiles %>%
filter(Data == "Buffalo" & variable == "canopy_mean"),
aes(x = hour, y = mean, colour = Data),
linewidth = 1) +
geom_path(data = hourly_summary_quantiles %>%
filter(!Data == "Buffalo" & variable == "canopy_mean"),
aes(x = hour, y = mean, colour = Data),
linewidth = 1) +
scale_fill_manual(values = colors) +
scale_colour_manual(values = colors) +
scale_x_discrete("Hour", breaks = seq(0,24,3)) +
scale_y_continuous("Mean value") +
ggtitle("Canopy cover") +
theme_classic()
hourly_path_canopy_plot
Code
hourly_path_slope_plot <- ggplot() +
geom_ribbon(data = hourly_summary_quantiles %>%
filter(Data == "Buffalo" & variable == "slope_mean"),
aes(x = hour, ymin = q25, ymax = q75, fill = Data),
alpha = ribbon_50_alpha) +
geom_ribbon(data = hourly_summary_quantiles %>%
filter(!Data == "Buffalo" & variable == "slope_mean"),
aes(x = hour, ymin = q25, ymax = q75, fill = Data),
alpha = ribbon_50_alpha) +
geom_path(data = summaries_hourly_all_long %>%
filter(Data == "Buffalo" & variable == "slope_mean"),
aes(x = hour, y = value, colour = Data, group = interaction(id, Data)),
alpha = buff_path_alpha,
linewidth = buff_path_linewidth) +
geom_path(data = hourly_summary_quantiles %>%
filter(!Data == "Buffalo" & variable == "slope_mean"),
aes(x = hour, y = q025, colour = Data),
linetype = "dashed",
alpha = path_95_alpha) +
geom_path(data = hourly_summary_quantiles %>%
filter(!Data == "Buffalo" & variable == "slope_mean"),
aes(x = hour, y = q975, colour = Data),
linetype = "dashed",
alpha = path_95_alpha) +
geom_path(data = hourly_summary_quantiles %>%
filter(Data == "Buffalo" & variable == "slope_mean"),
aes(x = hour, y = mean, colour = Data),
linewidth = 1) +
geom_path(data = hourly_summary_quantiles %>%
filter(!Data == "Buffalo" & variable == "slope_mean"),
aes(x = hour, y = mean, colour = Data),
linewidth = 1) +
scale_fill_manual(values = colors) +
scale_colour_manual(values = colors) +
scale_x_discrete("Hour", breaks = seq(0,24,3)) +
scale_y_continuous("Mean value") +
ggtitle("Slope") +
theme_classic()
hourly_path_slope_plot
Combining the hourly plots
Code
ggarrange(hourly_path_sl_plot +
theme(axis.title.x = element_blank(),
axis.text.x = element_blank()),
hourly_path_ndvi_plot +
theme(axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.title.y = element_blank()),
hourly_path_herby_plot,
hourly_path_canopy_plot +
theme(axis.title.y = element_blank()),
labels = c("A", "B", "C", "D"),
ncol = 2, nrow = 2,
legend = "bottom",
common.legend = TRUE)
References
Session info
R version 4.4.1 (2024-06-14 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 11 x64 (build 22631)
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
time zone: Europe/Zurich
tzcode source: internal
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] viridis_0.6.5 viridisLite_0.4.2 patchwork_1.2.0
[4] ggpubr_0.6.0 adehabitatHR_0.4.21 adehabitatLT_0.3.27
[7] CircStats_0.2-6 boot_1.3-30 MASS_7.3-60.2
[10] adehabitatMA_0.3.16 ade4_1.7-22 sp_2.1-4
[13] ks_1.14.2 beepr_2.0 tictoc_1.2.1
[16] terra_1.7-78 amt_0.2.2.0 lubridate_1.9.3
[19] forcats_1.0.0 stringr_1.5.1 dplyr_1.1.4
[22] purrr_1.0.2 readr_2.1.5 tidyr_1.3.1
[25] tibble_3.2.1 ggplot2_3.5.1 tidyverse_2.0.0
loaded via a namespace (and not attached):
[1] tidyselect_1.2.1 farver_2.1.2 fastmap_1.2.0 pracma_2.4.4
[5] digest_0.6.37 timechange_0.3.0 lifecycle_1.0.4 sf_1.0-16
[9] survival_3.6-4 magrittr_2.0.3 compiler_4.4.1 rlang_1.1.4
[13] tools_4.4.1 utf8_1.2.4 yaml_2.3.10 data.table_1.15.4
[17] knitr_1.48 ggsignif_0.6.4 labeling_0.4.3 htmlwidgets_1.6.4
[21] bit_4.0.5 mclust_6.1.1 classInt_0.4-10 abind_1.4-5
[25] KernSmooth_2.23-24 withr_3.0.1 grid_4.4.1 fansi_1.0.6
[29] e1071_1.7-14 colorspace_2.1-1 scales_1.3.0 cli_3.6.3
[33] mvtnorm_1.2-6 crayon_1.5.3 rmarkdown_2.28 generics_0.1.3
[37] rstudioapi_0.16.0 tzdb_0.4.0 DBI_1.2.3 proxy_0.4-27
[41] audio_0.1-11 splines_4.4.1 parallel_4.4.1 vctrs_0.6.5
[45] Matrix_1.7-0 jsonlite_1.8.8 carData_3.0-5 car_3.1-2
[49] hms_1.1.3 bit64_4.0.5 rstatix_0.7.2 units_0.8-5
[53] glue_1.7.0 codetools_0.2-20 cowplot_1.1.3 stringi_1.8.4
[57] gtable_0.3.5 munsell_0.5.1 pillar_1.9.0 htmltools_0.5.8.1
[61] R6_2.5.1 Rdpack_2.6.1 vroom_1.6.5 evaluate_0.24.0
[65] lattice_0.22-6 rbibutils_2.2.16 backports_1.5.0 broom_1.0.6
[69] class_7.3-22 Rcpp_1.0.13 checkmate_2.3.2 gridExtra_2.3
[73] xfun_0.47 pkgconfig_2.0.3