Manuscript Figure and Supplemental

Author

Steffi LaZerte

Published

September 23, 2024

This is the final figure and Supplemental material for the manuscript.

Setup

library(tidyverse)
library(patchwork)
library(ggspatial)
library(sf)
library(rnaturalearth)
library(ggthemes)
library(openxlsx)
library(gt)


source("XX_functions.R")

# Metrics
v <- read_csv("Data/Datasets/vultures_final.csv") |>
  # Round non-integer values of population counts
  mutate(across(c(contains("pop"), contains("raw")), round)) 

# Raw counts
raw <- read_csv("Data/Datasets/vultures_clean_2023.csv")

# Predicted GAM models
pred <- read_csv("Data/Datasets/vultures_gams_pred.csv")

# Checking problematic years
supp <- read_csv("Data/Datasets/table_supplemental.csv")

Main figure

Code
v <- v |> 
  mutate(date = as_date(p50_doy))

g1 <- ggplot(v, aes(x = year, y = date)) +
  theme_bw() +
  theme(axis.title.x = element_blank()) +
  geom_point() +
  geom_smooth(method = "lm", se = TRUE) +
  labs(y = "Date of peak migration", x = "")

g2 <- ggplot(v, aes(x = year, y = mig_raw_max)) +
  theme_bw() +
  geom_point() +
  stat_smooth(method = MASS::glm.nb) +
  labs(y = "Maximum observed count", x = "Year")

g <- g1 / g2 + plot_annotation(tag_levels = "A")
#ggsave("fig1_quick.png", dpi = 1000, width = 8, height = 7)

Big version

Code
g
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

Small version

Code
g
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

Map

Code
library(bcmaps)
library(ggrepel)
ne_download(type = "populated_places", scale = "large", load = FALSE)
[1] "ne_10m_populated_places"
Code
cities <- bcmaps::bc_cities(ask = FALSE) |> 
  select(NAME) |>
  filter(NAME %in% c("Victoria", "Vancouver", "Nanaimo"))
bc_cities was updated on 2024-07-10
Code
cities <- ne_load(file_name = "ne_10m_populated_places") |> 
  filter(NAME %in% c("Seattle", "Portland"), ADM1NAME %in% c("Oregon", "Washington")) |>
  select(NAME) |> 
  st_transform(st_crs(cities)) |>
  bind_rows(cities) |>
  rename(name = NAME)
Reading layer `ne_10m_populated_places' from data source 
  `/home/steffi/R_tmpdir/RtmpHNzeqe/ne_10m_populated_places.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 7342 features and 137 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -179.59 ymin: -90 xmax: 179.3833 ymax: 82.48332
Geodetic CRS:  WGS 84
Code
stn <- c(-123.55082035835214, 48.31773308537152) |>
  st_point() |>
  st_sfc(crs = 4326) |>
  st_sf(name = "Rocky Point\nObservatory") |>
  st_transform(st_crs(cities))

land <- data.frame(lon = c(-124.3, -124.3, -123.5),
                   lat = c(48.39, 48.75, 49.2),
                   name = c("Juan de Fuca Strait", "Vancouver Island", "Salish Sea")) |>
  st_as_sf(coords = c("lon", "lat"), crs = 4326)

area <- ne_states(country = c("Canada", "United States of America")) |>
  #st_crop(st_bbox(c(xmin = -140, xmax = -112, ymin = 41, ymax = 60))) |>
  #filter(name %in% c("British Columbia", "Alberta", "Washington", "Oregon", "Idaho", "Montana")) |>
  mutate(name = if_else(name %in% c("British Columbia", "Washington", "Oregon"), toupper(name), NA),
         name = str_replace(name, " ", "\n"),
         postal =  if_else(postal %in% c("BC", "WA", "OR", "AB"), postal, NA))

box <- st_polygon(list(rbind(c(1050000, 330000), 
                             c(1050000, 490000),
                             c(1230000, 490000),
                             c(1230000, 330000),
                             c(1050000, 330000)))) |>
  st_sfc(crs = 3005)

g0 <- ggplot(data = area, aes(label = name)) +
  theme_map() +
  theme(panel.border = element_rect(fill = NA), 
        plot.margin = unit(c(0,0,0,0), units = "mm"),
        panel.spacing = unit(0, units = "mm")) +
  scale_x_continuous(expand = c(0,0)) +
  scale_y_continuous(expand = c(0,0)) +
  geom_sf()
  
g <- g0 +
  geom_sf(data = stn, size = 3) +
  geom_sf_text(data = stn, lineheight = 0.85, nudge_y = -8000) +
  geom_sf(data = cities) +
  geom_sf_text(data = cities, hjust = 1.1, nudge_y = 5000) +
  geom_sf_text(data = land, angle = c(-24, 0, -45), colour = "grey60", size = c(6, 7, 7)) +
  annotation_scale(location = "bl") + 
  annotation_north_arrow(
    location = "br", height = unit(0.5, "cm"), width = unit(0.5, "cm"),
    style = north_arrow_orienteering(text_size = -Inf)) +
  coord_sf(crs = 3005, xlim = c(1050000, 1230000), ylim = c(330000, 490000))

g_inset <- g0 +
  theme(panel.background = element_rect(fill = "white")) +
  geom_sf_text(aes(label = postal), colour = "black", size = 3) +
  geom_sf(data = box, fill = NA, linewidth = 0.5, colour = "black", inherit.aes = FALSE) +
  coord_sf(crs = 3005, xlim = c(500000, 2000000), ylim = c(-200000, 1370000))

Big

Code
g + inset_element(g_inset, left = 0, top = 1, bottom = 0.6, right = 0.342, align_to = "full")
Warning: Removed 60 rows containing missing values or values outside the scale range
(`geom_text()`).

Small

Code
g + inset_element(g_inset, left = 0, top = 1, bottom = 0.6, right = 0.342, align_to = "full")
Warning: Removed 60 rows containing missing values or values outside the scale range
(`geom_text()`).

Supplemental Figure

Code
pred <- mutate(pred, ci95_upper = count + se * 1.96, ci95_lower = count - se * 1.96)
res <- filter(v, year == 1998) |>
  mutate(xmin = 204, 
         xmax = 240, y = 0.25 * mig_raw_max, 
         xmid = (240-204)/2 + 204,
         y = 0.25 * mig_raw_max,
         height = 0.03 * mig_raw_max, label = "Resident period")

g0 <- ggplot() +
  theme_void() +
  annotate(geom = "text", label = c("Migration\nstart", "Peak\nstart", "Peak", "Peak\nend", "Migration\nend"), 
           y = -0.25, x = c(0, 0.7, 1.3, 1.6, 2), size = 3, lineheight = 0.85) +
  annotate(geom = "text", label = c("5%", "25%", "50%", "75%", "95%"), y = 0.25, x = c(0, 0.7, 1.3, 1.6, 2), size = 3) +
  annotate(geom = "segment", y = 0, x = 0, xend = 2, arrow = arrow(angle = 90, ends = "both", length = unit(1, "mm"))) +
  annotate(geom = "segment", y = 0, x = 0.7, xend = 1.6, linewidth = 4) +
  annotate(geom = "segment", x = 1.29, xend = 1.30, y = 0, linewidth = 8) +
  ylim(c(-0.5, 0.5)) +
  xlim(c(-0.1, 2.1))

g1 <- ggplot(data = pred, mapping = aes(x = doy, y = count)) +
  theme_bw() +
  geom_ribbon(aes(ymin = ci95_lower, ymax = ci95_upper), fill = "grey50", alpha = 0.5) +
  geom_point(data = raw, na.rm = TRUE, size = 0.5) +
  geom_line() +
  geom_errorbarh(data = res, aes(xmin = xmin, xmax = xmax, y = y, height = height), 
                 colour = "grey70", inherit.aes = FALSE) +
  geom_text(data = res, aes(x = xmid, y = y, label = label), vjust = -0.5, inherit.aes = FALSE,
            size = 3, colour = "grey30") +
  geom_segment(data = v, aes(x = p50_doy, xend = p50_doy + 1, y = -(0.07 * mig_raw_max)), 
               linewidth = 4, inherit.aes = FALSE) +
  geom_errorbarh(data = v, aes(y = -(0.07 * mig_raw_max), 
                               xmin = mig_start_doy, 
                               xmax = mig_end_doy, 
                               height = 0.07 * mig_raw_max), inherit.aes = FALSE) +
  geom_segment(data = v, aes(y = -(0.07 * mig_raw_max), x = peak_start_doy, xend = peak_end_doy), linewidth = 2, inherit.aes = FALSE) +
  scale_x_continuous(name = "Date", limits = c(203, 295), 
                     labels = \(x) format(as_date(x), "%b %d"), 
                     n.breaks = 7) +
  labs(y = "Count") + 
  facet_wrap(~ year, scales = "free_y", ncol = 4)

Big

Code
g1 + inset_element(g0, left = 0.55, right = 0.95, bottom = 0, top = 0.1)

Small

Code
g1 + inset_element(g0, left = 0.55, right = 0.95, bottom = 0, top = 0.1)

Supplemental Table

Presented here for the record, but use exported XLSX version for submission.

Code
t <- supp |> 
  filter(!str_detect(term, "(I|i)ntercept"),
         !str_detect(model, "(max_doy)|(pop_max)")) |>
  mutate(across(where(is.numeric), \(x) round(x, 3))) |>
  mutate(p_value = if_else(p_value <= 0.05, 
                           paste0("**", format(p_value, nsmall = 3), "**"),
                           format(p_value, nsmall = 3)),
         value = paste0(format(estimate, nsmall = 3), " (P = ", p_value, ")"),
         value = str_trim(value)) |>
  select(Model = model, value, type) |>
  mutate(type = str_replace(type, "_", " "),
         type = str_to_title(type)) |>
  pivot_wider(names_from = "type", values_from = "value")

gt(t) |>
  gt_theme() |>
  fmt_markdown() |>
  fmt_number(columns = -1, decimals = 3)
Model Original Gaps Removed End Removed All Removed

mig_start_doy ~ year

0.176 (P = 0.017)

0.191 (P = 0.016)

0.166 (P = 0.037)

0.180 (P = 0.041)

peak_start_doy ~ year

0.155 (P = 0.005)

0.163 (P = 0.003)

0.155 (P = 0.014)

0.165 (P = 0.012)

p50_doy ~ year

0.147 (P = 0.009)

0.157 (P = 0.005)

0.151 (P = 0.024)

0.165 (P = 0.015)

peak_end_doy ~ year

0.132 (P = 0.016)

0.140 (P = 0.008)

0.128 (P = 0.043)

0.138 (P = 0.027)

mig_end_doy ~ year

0.093 (P = 0.263)

0.091 (P = 0.324)

0.077 (P = 0.402)

0.073 (P = 0.480)

mig_dur_days ~ year

-0.083 (P = 0.394)

-0.100 (P = 0.365)

-0.089 (P = 0.394)

-0.107 (P = 0.382)

peak_dur_days ~ year

-0.023 (P = 0.476)

-0.023 (P = 0.508)

-0.028 (P = 0.420)

-0.027 (P = 0.469)

mig_raw_max ~ year

0.041 (P = 0.000)

0.040 (P = 0.000)

0.034 (P = 0.002)

0.032 (P = 0.009)

Code
t <- supp |> 
  filter(!str_detect(term, "(I|i)ntercept")) |>
  mutate(across(where(is.numeric), \(x) round(x, 3))) |>
  select(model, type, estimate, p_value) |>
  pivot_wider(names_from = "type", values_from = c("estimate", "p_value"),
              names_glue = "{type}_{.value}") |>
  select(model, contains("original"), contains("gaps"), contains("end"), contains("all")) |>
  rename("Model" = "model",
         "Original" = "original_estimate", 
         "Removed Gaps" = "gaps_removed_estimate",
         "Removed Ends " = "end_removed_estimate", 
         "Removed All" = "all_removed_estimate")

wb <- createWorkbook()
addWorksheet(wb, "Supplemental Table 1")
writeData(wb, 1, t)
bold <- createStyle(textDecoration = "bold")
pval <- createStyle(halign = "left", numFmt = "(P = 0.000)")
est <- createStyle(halign = "right", numFmt = "0.000")
for(c in c(3, 5, 7, 9)) {
  conditionalFormatting(wb, 1, cols = c, rows = 1:15, rule = "<=0.05", style = bold)
  addStyle(wb, 1, cols = c, rows = 1:15, style = pval, stack = TRUE)
}
for(c in c(2, 4, 6, 8)) {
  addStyle(wb, 1, cols = c, rows = 1:15, style = est, stack = TRUE)
}
mergeCells(wb, 1, cols = 2:3, rows = 1)
mergeCells(wb, 1, cols = 4:5, rows = 1)
mergeCells(wb, 1, cols = 6:7, rows = 1)
mergeCells(wb, 1, cols = 8:9, rows = 1)
addStyle(wb, 1, cols = 1:9, rows = 1, 
         style = createStyle(halign = "center", textDecoration = "bold"),
         stack = TRUE)
setColWidths(wb, 1, cols = 1:9, widths = 10)
setColWidths(wb, 1, cols = 1, widths = 20)

setRowHeights(wb, 1, rows = nrow(t) + 3, heights = 100)
writeData(
  wb, 1, startCol = 1, startRow = nrow(t) + 3, 
  x = paste0(
    "Table. Results [Estimate of Year (P-value)] for different models when removing selective years. \n",
    "'Original' represents the results for the original model;\n",
    "'Removed Gaps' represents the results for these models when removing years with large gaps around the migration peaks (1998, 2011, 2013);\n",
    "'Removed Ends' represents the results for these models when removing years where it looks like sampling ended earlier than the end of migration (2001, 2006, 2023);\n",
    "'Removed All' represents the results for these models when removing all questionable years (1998, 2001, 2006, 2011, 2013, 2023)\n",
    "Bold P values are significant at P <= 0.05."
    ))

saveWorkbook(wb, "Data/table_supp.xlsx", overwrite = TRUE)
Back to top