KUD Habitat Comparisons
2026-04-22
library(dplyr)
library(trajr)
library(sf)
library(ggmap)
library(tidyverse)
library(rstanarm)
library(ggspatial)
library(gridExtra)
library(cowplot)
library(readr)
library(purrr)
library(lubridate)
library(zoo)
library(mapview)
library(emmeans)
library(geosphere)
library(slider)
library(scales)
library(ctmm)
mapviewOptions(fgb=F)habitats <- st_read('/Users/larryzheng/Desktop/Fulbright/Research/Drone Tracking/habitats_test4.geojson')## Reading layer `habitats_test4' from data source
## `/Users/larryzheng/Desktop/Fulbright/Research/Drone Tracking/habitats_test4.geojson'
## using driver `GeoJSON'
## Simple feature collection with 6 features and 1 field
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 103.0052 ymin: 5.812754 xmax: 103.0104 ymax: 5.81758
## Geodetic CRS: WGS 84mapview(habitats)Habitat overlaps for ARS
contour_level <- 0.5# Lets plot out what we have produced
cont_50_ars <- st_read('/Users/larryzheng/Desktop/Fulbright/Research/Drone Tracking/50_ars_contours.geojson')## Reading layer `50_contours' from data source
## `/Users/larryzheng/Desktop/Fulbright/Research/Drone Tracking/50_ars_contours.geojson'
## using driver `GeoJSON'
## Simple feature collection with 36 features and 2 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 103.0068 ymin: 5.812596 xmax: 103.0099 ymax: 5.814514
## Geodetic CRS: WGS 84mapview(cont_50_ars, zcol = "id", legend = F)Overlap the 50% KUD
overlaps_ars <-
st_intersection(cont_50_ars, habitats) %>%
mutate(area = as.numeric(st_area(.)))## Warning: attribute variables are assumed to be spatially constant throughout
## all geometriesoverlaps_ars %>%
mapview(zcol = "habitat")overlaps_ars %>%
st_drop_geometry() %>%
group_by(id, habitat) %>%
summarise(area = sum(area)) %>%
ggplot(aes(x = area, y = id, fill = habitat)) +
geom_col(position = "fill", col = "white", lwd = 0.1) +
labs(x = "Proportion of ARS activity space", fill = "Reef\nHabitat") +
scale_fill_brewer(palette = "Set1") +
scale_x_continuous(expand = c(0,0)) +
theme_bw()## `summarise()` has regrouped the output.
## ℹ Summaries were computed grouped by id and habitat.
## ℹ Output is grouped by id.
## ℹ Use `summarise(.groups = "drop_last")` to silence this message.
## ℹ Use `summarise(.by = c(id, habitat))` for per-operation grouping
## (`?dplyr::dplyr_by`) instead.
Habitat overlaps for Wandering
# Input wandering KUD files
cont_50_wandering <- st_read('/Users/larryzheng/Desktop/Fulbright/Research/Drone Tracking/50__wandering_contours.geojson')## Reading layer `50_contours' from data source
## `/Users/larryzheng/Desktop/Fulbright/Research/Drone Tracking/50__wandering_contours.geojson'
## using driver `GeoJSON'
## Simple feature collection with 36 features and 2 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 103.0055 ymin: 5.811821 xmax: 103.0102 ymax: 5.815312
## Geodetic CRS: WGS 84mapview(cont_50_wandering, zcol = "id", legend = F)Overlap the 50% KUD
overlaps_wandering <-
st_intersection(cont_50_wandering, habitats) %>%
mutate(area = as.numeric(st_area(.)))## Warning: attribute variables are assumed to be spatially constant throughout
## all geometriesoverlaps_wandering %>%
mapview(zcol = "habitat")overlaps_wandering %>%
st_drop_geometry() %>%
group_by(id, habitat) %>%
summarise(area = sum(area)) %>%
ggplot(aes(x = area, y = id, fill = habitat)) +
geom_col(position = "fill", col = "white", lwd = 0.1) +
labs(x = "Proportion of Wandering activity space", fill = "Reef\nHabitat") +
scale_fill_brewer(palette = "Set1") +
scale_x_continuous(expand = c(0,0)) +
theme_bw()## `summarise()` has regrouped the output.
## ℹ Summaries were computed grouped by id and habitat.
## ℹ Output is grouped by id.
## ℹ Use `summarise(.groups = "drop_last")` to silence this message.
## ℹ Use `summarise(.by = c(id, habitat))` for per-operation grouping
## (`?dplyr::dplyr_by`) instead.
overlaps_ars_test <- overlaps_ars %>%
group_by(id) %>%
mutate(total_area = sum(area),
prop = area / total_area) %>%
ungroup()
overlaps_wandering_test <- overlaps_wandering %>%
group_by(id) %>%
mutate(total_area = sum(area),
prop = area / total_area) %>%
ungroup()overlaps_ars_test$movement_type <- "ARS"
overlaps_wandering_test$movement_type <- "wandering"
overlaps_all <- bind_rows(overlaps_ars_test, overlaps_wandering_test)hist(overlaps_all$prop)
Fill missing habitats’ proportions with 0
overlaps_all_test <- overlaps_all %>%
complete(id, habitat, fill = list(area = 0, prop = 0))Proportions aren’t normal, so we must logit transform
epsilon <- 0.0001
overlaps_all_test <- overlaps_all_test %>%
mutate(prop_adj = pmin(pmax(prop, epsilon), 1 - epsilon),
prop_logit = log(prop_adj / (1 - prop_adj)))Clean and summarize dataframe
habitats <- c("lagoon", "reef_crest", "fore_reef", "sand_flat")
df_clean <- overlaps_all_test %>%
mutate(
habitat = tolower(habitat),
habitat = gsub(" ", "_", habitat),
prop = ifelse(is.na(prop), 0, prop)
) %>%
# keep only real habitat rows
filter(habitat %in% c("lagoon", "reef_crest", "fore_reef", "sand_flat")) %>%
group_by(id, movement_type, habitat, total_area) %>%
summarise(prop = sum(prop, na.rm = TRUE), .groups = "drop") %>%
pivot_wider(
names_from = habitat,
values_from = prop,
values_fill = 0
)
df_final <- df_clean %>%
mutate(
total = lagoon + reef_crest + fore_reef + sand_flat,
across(c(lagoon, reef_crest, fore_reef, sand_flat),
~ .x / total)
)Summary dataframe of means
df_final <- df_final %>%
filter(!is.na(movement_type))
df_final %>%
group_by(movement_type) %>%
summarise(across(all_of(habitats), mean))## # A tibble: 2 × 5
## movement_type lagoon reef_crest fore_reef sand_flat
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 ARS 0.887 0.0193 0.0159 0.0777
## 2 wandering 0.575 0.0418 0.139 0.245Summary dataframe of both means and sd
df_summary <- df_final %>%
group_by(movement_type) %>%
summarise(
across(
all_of(habitats),
list(
mean = ~ mean(.x, na.rm = TRUE),
sd = ~ sd(.x, na.rm = TRUE)
),
.names = "{.col}_{.fn}"
)
)
df_summary## # A tibble: 2 × 9
## movement_type lagoon_mean lagoon_sd reef_crest_mean reef_crest_sd
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 ARS 0.887 0.232 0.0193 0.0495
## 2 wandering 0.575 0.420 0.0418 0.0527
## # ℹ 4 more variables: fore_reef_mean <dbl>, fore_reef_sd <dbl>,
## # sand_flat_mean <dbl>, sand_flat_sd <dbl>PERMANOVA
library(vegan)## Loading required package: permuteadonis2(df_final[, habitats] ~ movement_type,
data = df_final,
method = "bray",
permutations = 9999)## Permutation test for adonis under reduced model
## Permutation: free
## Number of permutations: 9999
##
## adonis2(formula = df_final[, habitats] ~ movement_type, data = df_final, permutations = 9999, method = "bray")
## Df SumOfSqs R2 F Pr(>F)
## Model 1 1.5910 0.15156 12.504 8e-04 ***
## Residual 70 8.9069 0.84844
## Total 71 10.4980 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Bray-Curtis dissimilarity
d <- vegdist(df_final[, habitats], method = "bray")
# PERMANOVA model
adonis_model <- adonis2(d ~ movement_type, data = df_final)
# check variable contributions via envfit
fit <- envfit(prcomp(df_final[, habitats]), df_final$movement_type)
fit##
## ***FACTORS:
##
## Centroids:
## PC1 PC2
## df_final.movement_typeARS -0.1856 -0.0291
## df_final.movement_typewandering 0.1856 0.0291
##
## Goodness of fit:
## r2 Pr(>r)
## df_final.movement_type 0.1519 0.001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999SIMPER Test (finds that lagoon (reef flat) is biggest contributor to difference)
simper(df_final[, habitats], df_final$movement_type)## cumulative contributions of most influential species:
##
## $ARS_wandering
## lagoon sand_flat
## 0.4813165 0.7825026PERMANOVA controlling for ID
adonis2(df_final[, habitats] ~ movement_type,
data = df_final,
method = "bray",
permutations = 9999,
strata = df_final$id)## Permutation test for adonis under reduced model
## Blocks: strata
## Permutation: free
## Number of permutations: 9999
##
## adonis2(formula = df_final[, habitats] ~ movement_type, data = df_final, permutations = 9999, method = "bray", strata = df_final$id)
## Df SumOfSqs R2 F Pr(>F)
## Model 1 1.5910 0.15156 12.504 0.2017
## Residual 70 8.9069 0.84844
## Total 71 10.4980 1.00000GLMM
library(glmmTMB)## Warning in check_dep_version(dep_pkg = "TMB"): package version mismatch:
## glmmTMB was built with TMB package version 1.9.19
## Current TMB package version is 1.9.21
## Please re-install glmmTMB from source or restore original 'TMB' package (see '?reinstalling' for more information)epsilon <- 0.001
overlaps_all_test2 <- overlaps_all_test %>%
mutate(prop_adj = pmin(pmax(prop, epsilon), 1 - epsilon))
glmm_model <- glmmTMB(prop_adj ~ movement_type * habitat + (1 | id),
family = beta_family(), data = overlaps_all_test2)emmeans(glmm_model, pairwise ~ movement_type | habitat)## $emmeans
## habitat = fore reef:
## movement_type emmean SE df asymp.LCL asymp.UCL
## ARS -1.186 0.462 Inf -2.091 -0.2805
## wandering -1.280 0.214 Inf -1.699 -0.8609
##
## habitat = lagoon:
## movement_type emmean SE df asymp.LCL asymp.UCL
## ARS 2.048 0.203 Inf 1.651 2.4456
## wandering 1.229 0.207 Inf 0.823 1.6360
##
## habitat = reef crest:
## movement_type emmean SE df asymp.LCL asymp.UCL
## ARS -1.424 0.399 Inf -2.207 -0.6416
## wandering -1.599 0.212 Inf -2.015 -1.1827
##
## habitat = sand flat:
## movement_type emmean SE df asymp.LCL asymp.UCL
## ARS -0.759 0.380 Inf -1.504 -0.0143
## wandering -0.316 0.249 Inf -0.804 0.1721
##
## Results are given on the logit (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## habitat = fore reef:
## contrast estimate SE df z.ratio p.value
## ARS - wandering 0.0945 0.501 Inf 0.188 0.8505
##
## habitat = lagoon:
## contrast estimate SE df z.ratio p.value
## ARS - wandering 0.8186 0.271 Inf 3.022 0.0025
##
## habitat = reef crest:
## contrast estimate SE df z.ratio p.value
## ARS - wandering 0.1743 0.441 Inf 0.396 0.6924
##
## habitat = sand flat:
## contrast estimate SE df z.ratio p.value
## ARS - wandering -0.4427 0.453 Inf -0.978 0.3282
##
## Results are given on the log odds ratio (not the response) scale.emmeans(glmm_model, ~ movement_type * habitat, type = "response")## movement_type habitat response SE df asymp.LCL asymp.UCL
## ARS fore reef 0.234 0.0828 Inf 0.1100 0.430
## wandering fore reef 0.218 0.0364 Inf 0.1545 0.297
## ARS lagoon 0.886 0.0205 Inf 0.8390 0.920
## wandering lagoon 0.774 0.0363 Inf 0.6948 0.837
## ARS reef crest 0.194 0.0624 Inf 0.0991 0.345
## wandering reef crest 0.168 0.0297 Inf 0.1177 0.235
## ARS sand flat 0.319 0.0825 Inf 0.1819 0.496
## wandering sand flat 0.422 0.0608 Inf 0.3091 0.543
##
## Confidence level used: 0.95
## Intervals are back-transformed from the logit scaleGraph results of GLMM
Inputted values from above manually
# Create data
df <- data.frame(
Habitat = rep(c("Fore reef", "Lagoon", "Reef crest", "Sand flat"), each = 2),
Movement = rep(c("ARS", "Wandering"), times = 4),
Mean = c(0.234, 0.218,
0.886, 0.774,
0.194, 0.168,
0.319, 0.422),
CI_low = c(0.110, 0.155,
0.839, 0.695,
0.099, 0.118,
0.182, 0.309),
CI_high = c(0.430, 0.297,
0.920, 0.837,
0.345, 0.235,
0.496, 0.543),
p_value = rep(c(0.8505, 0.0025, 0.6924, 0.3282), each = 2)
)
# Flag lagoon for highlighting (only lagoon has statistically significant p value)
df <- df %>%
mutate(Highlight = ifelse(Habitat == "Lagoon", "Lagoon", "Other"))ggplot(df, aes(x = Habitat, y = Mean, color = Movement)) +
# Points
geom_point(aes(alpha = Highlight),
size = 3.5,
stroke = 1,
position = position_dodge(width = 0.5)) +
# Error bars
geom_errorbar(aes(ymin = CI_low, ymax = CI_high, alpha = Highlight),
width = 0.2,
linewidth = 1,
position = position_dodge(width = 0.5)) +
# Custom colors
scale_color_manual(
name = "Movement Type",
values = c("ARS" = "#d73027", "Wandering" = "#4575b4")
) +
# Highlight lagoon via alpha
scale_alpha_manual(values = c("Lagoon" = 1, "Other" = 0.5), guide = "none") +
labs(y = "Mean Proportion", x = "Habitat") +
theme_minimal() +
theme(
axis.line.x = element_line(color = "black", linewidth = 0.6),
axis.line.y = element_line(color = "black", linewidth = 0.6),
panel.grid.major.x = element_blank(),
panel.grid.minor = element_blank(),
panel.grid.major.y = element_line(color = "gray80", linewidth = 0.3),
axis.title.x = element_text(size = 17),
axis.title.y = element_text(size = 17),
axis.text.x = element_text(size = 16),
axis.text.y = element_text(size = 16),
legend.title = element_text(size = 16),
legend.text = element_text(size = 16),
legend.position = c(0.85, 0.85),
legend.background = element_rect(fill = "white", color = NA)
)
GLM
library(betareg)
glm_model <- betareg(prop_adj ~ movement_type * habitat,
data = overlaps_all_test2)
emmeans(glm_model, pairwise ~ movement_type | habitat)## $emmeans
## habitat = fore reef:
## movement_type emmean SE df asymp.LCL asymp.UCL
## ARS 0.234 0.0828 Inf 0.0718 0.396
## wandering 0.218 0.0364 Inf 0.1462 0.289
##
## habitat = lagoon:
## movement_type emmean SE df asymp.LCL asymp.UCL
## ARS 0.886 0.0205 Inf 0.8455 0.926
## wandering 0.774 0.0363 Inf 0.7025 0.845
##
## habitat = reef crest:
## movement_type emmean SE df asymp.LCL asymp.UCL
## ARS 0.194 0.0624 Inf 0.0716 0.316
## wandering 0.168 0.0297 Inf 0.1100 0.226
##
## habitat = sand flat:
## movement_type emmean SE df asymp.LCL asymp.UCL
## ARS 0.319 0.0825 Inf 0.1572 0.481
## wandering 0.422 0.0608 Inf 0.3025 0.541
##
## Confidence level used: 0.95
##
## $contrasts
## habitat = fore reef:
## contrast estimate SE df z.ratio p.value
## ARS - wandering 0.0165 0.0891 Inf 0.185 0.8530
##
## habitat = lagoon:
## contrast estimate SE df z.ratio p.value
## ARS - wandering 0.1120 0.0394 Inf 2.847 0.0044
##
## habitat = reef crest:
## contrast estimate SE df z.ratio p.value
## ARS - wandering 0.0258 0.0675 Inf 0.383 0.7021
##
## habitat = sand flat:
## contrast estimate SE df z.ratio p.value
## ARS - wandering -0.1027 0.1020 Inf -1.006 0.3144