library(tidyverse)
library(patchwork)
library(ggspatial)
library(sf)
library(rnaturalearth)
library(ggthemes)
library(openxlsx)
library(gt)
source("XX_functions.R")
# Metrics
<- read_csv("Data/Datasets/vultures_final.csv") |>
v # Round non-integer values of population counts
mutate(across(c(contains("pop"), contains("raw")), round))
# Raw counts
<- read_csv("Data/Datasets/vultures_clean_2023.csv")
raw
# Predicted GAM models
<- read_csv("Data/Datasets/vultures_gams_pred.csv")
pred
# Checking problematic years
<- read_csv("Data/Datasets/table_supplemental.csv") supp
Manuscript Figure and Supplemental
This is the final figure and Supplemental material for the manuscript.
Setup
Main figure
Code
<- v |>
v mutate(date = as_date(p50_doy) - days(1))
<- ggplot(v, aes(x = year, y = date)) +
g1 theme_bw() +
theme(axis.title.x = element_blank()) +
geom_point() +
geom_smooth(method = "lm", se = TRUE) +
labs(y = "Date of peak migration", x = "")
<- ggplot(v, aes(x = year, y = mig_raw_max)) +
g2 theme_bw() +
geom_point() +
stat_smooth(method = MASS::glm.nb) +
labs(y = "Annual maximum\nDaily Estimated Total", x = "Year")
<- g1 / g2 + plot_annotation(tag_levels = "A")
g #ggsave("fig1_quick.png", dpi = 1000, width = 8, height = 7)
<- g1 / (g2 + labs(y = "Annual maximum DET")) + plot_annotation(tag_levels = "A") gg
Double checking growth and increases
Code
<- MASS::glm.nb(mig_raw_max ~ year, data = v)
m
<- select(v, year) |>
d mutate(y = predict(m, v),
yexp = exp(y))
# Compound growth (avg growth over the period)
1251 / 451)^ (1/25) - 1 (
[1] 0.04165339
Code
# Compound interest (amount at the end)
451 * (1 + 0.041647)^(25)
[1] 1250.808
Code
# Factor of increase
1251/451
[1] 2.773836
Big version
With two y-axis options for panel B
Code
g
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
Code
gg
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
Small version
With two y-axis options for panel B
Code
g
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
Code
gg
`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
<- bcmaps::bc_cities(ask = FALSE) |>
cities select(NAME) |>
filter(NAME %in% c("Victoria", "Vancouver", "Nanaimo"))
bc_cities was updated on 2024-07-10
Code
<- ne_load(file_name = "ne_10m_populated_places") |>
cities 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/RtmpTPJOlj/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
<- c(-123.55082035835214, 48.31773308537152) |>
stn st_point() |>
st_sfc(crs = 4326) |>
st_sf(name = "Rocky Point") |>
st_transform(st_crs(cities))
<- data.frame(lon = c(-124.3, -124.3, -123.5),
land 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)
<- ne_states(country = c("Canada", "United States of America")) |>
area #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))
<- st_polygon(list(rbind(c(1050000, 330000),
box c(1050000, 490000),
c(1230000, 490000),
c(1230000, 330000),
c(1050000, 330000)))) |>
st_sfc(crs = 3005)
<- ggplot(data = area, aes(label = name)) +
g0 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()
<- g0 +
g 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))
<- g0 +
g_inset 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
+ inset_element(g_inset, left = 0, top = 1, bottom = 0.6, right = 0.342, align_to = "full") g
Warning: Removed 60 rows containing missing values or values outside the scale range
(`geom_text()`).
Small
Code
+ inset_element(g_inset, left = 0, top = 1, bottom = 0.6, right = 0.342, align_to = "full") g
Warning: Removed 60 rows containing missing values or values outside the scale range
(`geom_text()`).
Supplemental Figure
Code
<- mutate(pred, ci95_upper = count + se * 1.96, ci95_lower = count - se * 1.96)
pred <- filter(v, year == 1998) |>
res 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")
<- ggplot() +
g0 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))
<- ggplot(data = pred, mapping = aes(x = doy, y = count)) +
g1 theme_bw() +
# GAM
geom_ribbon(aes(ymin = ci95_lower, ymax = ci95_upper), fill = "grey50", alpha = 0.5) +
geom_line() +
# Raw points
geom_point(data = raw, na.rm = TRUE, size = 0.5) +
# Metrics
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 - 0.5, xend = p50_doy + 0.5,
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) - days(1), "%b %d"),
n.breaks = 7) +
labs(y = "Daily Estimated Total") +
facet_wrap(~ year, scales = "free_y", ncol = 4)
Big
Code
+ inset_element(g0, left = 0.55, right = 0.95, bottom = 0, top = 0.1) g1
Small
Code
+ inset_element(g0, left = 0.55, right = 0.95, bottom = 0, top = 0.1) g1
Supplemental Table
Presented here for the record, but use exported XLSX version for submission.
Code
<- supp |>
t 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
<- supp |>
t 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")
<- createWorkbook()
wb addWorksheet(wb, "Supplemental Table 1")
writeData(wb, 1, t)
<- createStyle(textDecoration = "bold")
bold <- createStyle(halign = "left", numFmt = "(P = 0.000)")
pval <- createStyle(halign = "right", numFmt = "0.000")
est 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(
1, startCol = 1, startRow = nrow(t) + 3,
wb, 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)