---
title: "Evaluating the effects of tracking devices on survival, breeding success, behaviour, and condition of a small, partially migratory shorebird"
subtitle: "Supplementary Material A: Journal of Avian Biology 10.1002/jav.03490"
date: "`r format(Sys.time(), '%d %B, %Y')`"
author:
- name: Luke Eberhart-Hertel
orcid: 0000-0001-7311-6088
email: luke.eberhart@bi.mpg.de
url: https://www.bi.mpg.de/person/115852/2867
affiliations:
- ref: bk
- name: Emma Williams
affiliations:
- ref: ew
- name: Ailsa McGilvary-Howard
affiliations:
- ref: ah
- name: Ted Howard
affiliations:
- ref: ah
- name: Tony Habraken
affiliations:
- ref: th
- name: Colin O`Donnell
affiliations:
- ref: ew
- name: Clemens Küpper
affiliations:
- ref: ck
- name: Bart Kempenaers
affiliations:
- ref: bk
affiliations:
- id: bk
number: 1
name: Department of Ornithology, Max Planck Institute for Biological Intelligence, Eberhard-Gwinner Straße 7, 82319 Seewiesen, Germany
- id: ah
number: 2
name: Kaikōura Banded Dotterel Project, 1 Maui Street, Kaikōura, New Zealand
- id: th
number: 2
name: Port Waikato Banded Dotterel Project, Port Waikato, New Zealand
- id: ew
number: 3
name: Fauna Science Team, Department of Conservation Christchurch Office, Christchurch Mail Centre, Private Bag 4715, Christchurch, 8140, New Zealand
- id: ck
number: 1
name: Behavioural Genetics and Evolutionary Ecology, Max Planck Institute for Biological Intelligence, Eberhard-Gwinner Straße 5, 82319 Seewiesen, Germany
format:
html:
toc: true
code-fold: false
code-tools: true
self-contained: true
highlight-style: github
theme: Cosmo
execute:
warning: false
cache: true
freeze: auto
editor_options:
chunk_output_type: console
---
```{r, include=FALSE}
knitr::opts_chunk$set(cache = TRUE)
```
## Prerequisites
### R packages
```{r, message=FALSE, results='hide', warning=FALSE, cache=FALSE}
# Set a CRAN mirror before installing packages
options(repos = c(CRAN = "https://cloud.r-project.org"))
# a vector of all the packages needed in the project
packages_required_in_project <-
c("apis", "BaSTA", "brms", "broom.mixed", "colorspace", "corrplot", "data.table",
"effects", "elevatr", "geosphere", "ggblend", "ggeffects", "ggmap",
"ggnewscale", #"ggsn",
"giscoR", "glue", "googlesheets4", "gt", "gtsummary",
"here", "hms", "lme4", "lubridate", "mapview", "mgcv", "multcomp", "MuMIn", "nnet",
"partR2", "patchwork", "RColorBrewer", "readxl", "remotes", "RMark", "rnaturalearth",
"rptR", "scales", "sf", "showtext", "smatr", "sp", "stringr", "tidybayes",
"tidyterra", "tidyverse", "adehabitatLT", "magrittr", "MASS", "leaflet", "shiny",
"leaflet.extras", "ggrepel", "ggpmisc", "grid", "knitr", "kableExtra")
# of the required packages, check if some need to be installed
new.packages <-
packages_required_in_project[!(packages_required_in_project %in%
installed.packages()[,"Package"])]
# install all packages that are not locally available
if(length(new.packages)) install.packages(new.packages)
# load all the packages into the current R session
lapply(packages_required_in_project, require, character.only = TRUE)
```
```{r, include=FALSE}
# set the home directory to where the project is locally based (i.e., to find
# the relevant datasets to import, etc.
here::set_here()
library(DBI)
library(dbo)
# Use gs4_auth() to log in and authenticate with your Google account.
# You can specify your email address to ensure the correct account is used.
# Once authenticated, you won't need to select the account interactively when using read_sheet()
gs4_auth(email = "luke.eberhart@gmail.com")
```
```{r, message=FALSE, results='hide', warning=FALSE, include=FALSE}
### Plotting themes
# - The following plotting themes, colors, and typefaces are used throughout the project:
# Find fonts from computer that you want. Use regular expressions to do this
# For example, load all fonts that are 'verdana' or 'Verdana'
extrafont::font_import(pattern = "[V/v]erdana", prompt = FALSE)
# extrafont::font_import(pattern = "[L/l]ato", prompt = FALSE)
# check which fonts were loaded
extrafont::fonts()
extrafont::fonttable()
extrafont::loadfonts() # load these into R
# define the plotting theme to be used in subsequent ggplots
luke_theme <-
theme_bw() +
theme(
text = element_text(family = "Verdana"),
legend.title = element_text(size = 10),
legend.text = element_text(size = 8),
axis.title.x = element_text(size = 10),
axis.text.x = element_text(size = 8),
axis.title.y = element_text(size = 10),
axis.text.y = element_text(size = 8),
strip.text = element_text(size = 10),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.ticks = element_line(linewidth = 0.3, colour = "grey40"),
axis.ticks.length = unit(0.2, "cm"),
panel.border = element_rect(linetype = "solid", colour = "grey"),
legend.position.inside = c(0.1, 0.9)
)
```
### Custom functions
```{r}
# function that does the opposite of "%in%"
` %!in% ` = Negate (` %in% ` )
# scaled mass index calculation from Peig and Green (2009)
# https://doi.org/10.1111/j.1365-2435.2010.01751.x
scaledMassIndex <-
function (x, y, x.0 = mean (x)) {
logM.ols <- lm (log (y) ~ log (x))
logM.rob <- rlm (log (y) ~ log (x), method = "M" )
b.msa.ols <- coef (sma (log (y) ~ log (x)))[2 ]
b.msa.rob <- coef (sma (log (y) ~ log (x), robust = T))[2 ]
SMI.ols <- y * (x.0 / x) ^ b.msa.ols
SMI.rob <- y * (x.0 / x) ^ b.msa.rob
res <- data.frame (SMI.ols, SMI.rob, x, y)
pred.DT <-
data.table (x = seq (min (x), max (x), length = 100 )) %>%
.[, y.ols : = predict (logM.ols, newdata = .) %>% exp] %>%
.[, y.rob : = predict (logM.rob, newdata = .) %>% exp]
attr (res, "b.msa" ) <- c (ols = b.msa.ols, rob = b.msa.rob)
return (res)
}
```
```{r, echo=FALSE}
# convert sf geometry to columns
sfc_as_cols <- function(x, geometry, names = c("x","y")) {
if (missing(geometry)) {
geometry <- sf::st_geometry(x)
} else {
geometry <- rlang::eval_tidy(enquo(geometry), x)
}
stopifnot(inherits(x,"sf") && inherits(geometry,"sfc_POINT"))
ret <- sf::st_coordinates(geometry)
ret <- tibble::as_tibble(ret)
stopifnot(length(names) == ncol(ret))
x <- x[ , !names(x) %in% names]
ret <- setNames(ret,names)
dplyr::bind_cols(x,ret)
}
```
# Body Condition Dynamics
### specify attachment weights
```{r}
harness = 0.29
colour_band = 0.08
metal_band = 0.16
PTT = 1.8
PTT2 = 2.0
GPS = 1.2
```
## PCA of structural traits
```{r}
body_traits <-
read.csv (here ("data/body_traits_season.csv" ))
```
#### check most relevent structural traits with PCA
```{r}
body_traits_pca <-
body_traits %>%
group_by (ring) %>%
arrange (measure) %>%
slice (1 ) %>%
ungroup () %>%
dplyr:: select (culmen, tarsus, wing) %>%
na.omit () %>%
princomp ()
# check the PCA results
summary (body_traits_pca)
biplot (body_traits_pca, cex = 0.7 )
loadings (body_traits_pca)
# bind PC scores to the original dataframe
body_traits_pca_scores <-
body_traits %>%
group_by (ring) %>%
arrange (measure) %>%
slice (1 ) %>%
bind_cols (., body_traits_pca$ scores[, 1 ], body_traits_pca$ scores[, 2 ], body_traits_pca$ scores[, 3 ]) %>%
rename (structure_pc1 = ` ...13 ` ,
structure_pc2 = ` ...14 ` ,
structure_pc3 = ` ...15 ` ) %>%
na.omit () %>%
ungroup ()
```
#### check correlations of traits with weight to determine best trait to use for SMI
```{r, echo=FALSE}
body_traits_pca_scores %>%
mutate(log_weight = log(weight),
log_culmen = log(culmen),
log_tarsus = log(tarsus),
log_wing = log(wing)) %>%
dplyr::select(log_weight, log_culmen, log_tarsus, log_wing, structure_pc1, structure_pc2, structure_pc3) %>%
cor() %>%
corrplot(type = "upper", method = "number", tl.srt = 45)
```
#### visualize the relationship between the principle components and body mass
```{r, echo=FALSE}
ggplot(body_traits_pca_scores,
aes(x = structure_pc1, y = weight)) +
geom_point() +
labs(x = "PC1", y = "Body Mass (g)") +
theme_minimal()
ggplot(body_traits_pca_scores,
aes(x = structure_pc2, y = weight)) +
geom_point() +
labs(x = "PC2", y = "Body Mass (g)") +
theme_minimal()
ggplot(body_traits_pca_scores,
aes(x = structure_pc3, y = weight)) +
geom_point() +
labs(x = "PC3", y = "Body Mass (g)") +
theme_minimal()
```
#### visualize the relationship between the wing length and body mass
```{r, echo=FALSE}
ggplot(body_traits_pca_scores,
aes(x = wing, y = weight)) +
geom_point() +
labs(x = "wing length (mm)", y = "Body Mass (g)") +
theme_minimal()
```
#### visualize the relationship between the tarsus length and body mass
```{r, echo=FALSE}
ggplot(body_traits_pca_scores,
aes(x = wing, y = tarsus)) +
geom_point() +
labs(x = "tarsus length (mm)", y = "Body Mass (g)") +
theme_minimal()
```
## Scaled Mass Index wrangle
average all repeated measures of structural traits for each bird (i.e., variation within individuals for these traits is assumed to be observer/instrument error)
```{r}
body_traits_avg <-
body_traits %>%
group_by (ring) %>%
mutate (culmen = mean (culmen, na.rm = TRUE ),
tarsus = mean (tarsus, na.rm = TRUE ),
wing = mean (wing, na.rm = TRUE ))
```
#### determine which structural trait to use for the scaled mass index calculation
see text in Peig and Green (2009) after Eq. 2: "...We recommend the use of the single 'L' variable which has the strongest correlation with M on a log-log scale, since this is likely to be the L that best explains that fraction of mass associated with structural size"
Conclude to use wing for SMI transformation because it has the largest r and based on the PCA investigation above, wing did just as good as PC1 when explaining variation in weight
```{r}
cor.test (log (body_traits_avg$ tarsus), log (body_traits_avg$ weight)) # r = 0.130
cor.test (log (body_traits_avg$ wing), log (body_traits_avg$ weight)) # r = 0.226
cor.test (log (body_traits_avg$ culmen), log (body_traits_avg$ weight)) # r = 0.090
```
Use scaledMassIndex() function [ defined in chunk above ] to derive the SMI for each body mass measure (see Eq 2 of Peig and Green (2009))
```{r}
smi_ <-
scaledMassIndex (x = body_traits_avg %>%
ungroup () %>%
dplyr:: select (weight, wing) %>%
na.omit () %>%
pull (wing),
y = body_traits_avg %>%
ungroup () %>%
dplyr:: select (weight, wing) %>%
na.omit () %>%
pull (weight))
body_traits_smi <-
body_traits_avg %>%
filter (! is.na (weight)) %>%
bind_cols (., smi_) %>%
dplyr:: select (ring, tag_type, tag_type_, treatment_group, measure, SMI.ols, weight, day_of_year, date_diff) %>%
rename (smi_wing = SMI.ols)
```
inspect relationship between raw body mass and scaled mass index (lines connect repeated measures)
```{r, echo=FALSE}
ggplot(body_traits_smi,
aes(x = smi_wing, y = weight)) +
geom_point(alpha = 0.2) +
geom_line(aes(group = ring, color = tag_type), linewidth = 1) +
labs(x = "scaled mass index (g-hat)", y = "raw body mass (g)") +
theme_minimal() +
scale_color_brewer(palette = "Dark2", name = "tag type")
```
## Scaled Mass Index sequential measures model
#### sample size summary
```{r, echo=FALSE}
body_traits_smi %>%
group_by(tag_type) %>%
summarise(n = n_distinct(ring))
```
#### summary of attachment mass
```{r, echo=FALSE}
# Function to replace duplicate values with blanks
remove_redundant <- function(x) {
dplyr::if_else(duplicated(x), "", x)
}
body_traits_smi %>%
group_by(ring) %>%
arrange(measure) %>%
slice(1) %>%
ungroup() %>%
mutate(prop_mass_tag = ifelse(tag_type_ == "PTT2", PTT2 / weight,
ifelse(tag_type_ == "PTT", PTT / weight,
ifelse(tag_type_ == "GPS", GPS / weight,
ifelse(tag_type_ == "control", NA, "XXX")))) %>% as.numeric(),
prop_mass_all = ifelse(tag_type_ == "PTT2", (PTT2 + sum(4 * colour_band, metal_band)) / weight,
ifelse(tag_type_ == "PTT", (PTT + sum(4 * colour_band, metal_band)) / weight,
ifelse(tag_type_ == "GPS", (GPS + sum(4 * colour_band, metal_band)) / weight,
ifelse(tag_type_ == "control", sum(4 * colour_band, metal_band) / weight, "XXX")))) %>% as.numeric(),
prop_SMI_tag = ifelse(tag_type_ == "PTT2", PTT2 / smi_wing,
ifelse(tag_type_ == "PTT", PTT / smi_wing,
ifelse(tag_type_ == "GPS", GPS / smi_wing,
ifelse(tag_type_ == "control", NA, "XXX")))) %>% as.numeric(),
prop_SMI_all = ifelse(tag_type_ == "PTT2", (PTT2 + sum(4 * colour_band, metal_band)) / smi_wing,
ifelse(tag_type_ == "PTT", (PTT + sum(4 * colour_band, metal_band)) / smi_wing,
ifelse(tag_type_ == "GPS", (GPS + sum(4 * colour_band, metal_band)) / smi_wing,
ifelse(tag_type_ == "control", sum(4 * colour_band, metal_band) / smi_wing, "XXX")))) %>% as.numeric()) %>%
ungroup() %>%
group_by(tag_type) %>%
summarise(mean_tag_SMI = mean(prop_SMI_tag, na.rm = TRUE) * 100,
median_tag_SMI = median(prop_SMI_tag, na.rm = TRUE) * 100,
min_tag_SMI = min(prop_SMI_tag, na.rm = TRUE) * 100,
max_tag_SMI = max(prop_SMI_tag, na.rm = TRUE) * 100,
sd_tag_SMI = sd(prop_SMI_tag, na.rm = TRUE) * 100,
mean_all_SMI = mean(prop_SMI_all, na.rm = TRUE) * 100,
median_all_SMI = median(prop_SMI_all, na.rm = TRUE) * 100,
min_all_SMI = min(prop_SMI_all, na.rm = TRUE) * 100,
max_all_SMI = max(prop_SMI_all, na.rm = TRUE) * 100,
sd_all_SMI = sd(prop_SMI_all, na.rm = TRUE) * 100,
mean_tag_mass = mean(prop_mass_tag, na.rm = TRUE) * 100,
median_tag_mass = median(prop_mass_tag, na.rm = TRUE) * 100,
min_tag_mass = min(prop_mass_tag, na.rm = TRUE) * 100,
max_tag_mass = max(prop_mass_tag, na.rm = TRUE) * 100,
sd_tag_mass = sd(prop_mass_tag, na.rm = TRUE) * 100,
mean_all_mass = mean(prop_mass_all, na.rm = TRUE) * 100,
median_all_mass = median(prop_mass_all, na.rm = TRUE) * 100,
min_all_mass = min(prop_mass_all, na.rm = TRUE) * 100,
max_all_mass = max(prop_mass_all, na.rm = TRUE) * 100,
sd_all_mass = sd(prop_mass_all, na.rm = TRUE) * 100
) %>%
t() %>% as.data.frame() %>% rename(GPS = V1,
PTT = V2,
# PTT2 = V3,
control = V3) %>% slice(-1) %>%
rownames_to_column(var = "statistic") %>%
mutate(GPS = ifelse(!is.infinite(GPS) | !is.na(GPS), floor(as.numeric(GPS) * 10) / 10, NA),
PTT = ifelse(!is.infinite(PTT) | !is.na(PTT), floor(as.numeric(PTT) * 10) / 10, NA),
control = ifelse(!is.infinite(control) | !is.na(control), floor(as.numeric(control) * 10) / 10, NA)) %>%
mutate(unit = ifelse(str_detect(statistic, "SMI"), "scaled mass index", "raw body mass"),
attachment = ifelse(str_detect(statistic, "tag"), "only tracking device", "all attachments"),
statistic = ifelse(str_detect(statistic, "mean"), "mean",
ifelse(str_detect(statistic, "median"), "median",
ifelse(str_detect(statistic, "min"), "min",
ifelse(str_detect(statistic, "max"), "max", "sd"))))) %>%
mutate(control = ifelse(is.infinite(control), NA, control)) %>%
group_by(unit) %>%
mutate(across(c(attachment), remove_redundant)) %>%
dplyr::select(unit, attachment, statistic, GPS, PTT, control) %>%
mutate(unit = factor(unit, levels = c("raw body mass", "scaled mass index"))) %>%
arrange(unit) %>%
gt(rowname_col = "attachment",
groupname_col = "unit") %>%
cols_label(statistic = "",
GPS = "1.2 g GPS (%)",
PTT = "1.8-2.0 g PTT (%)",
control = "only leg bands (%)") %>%
cols_align(
align = "center",
columns = c(GPS, PTT, control)
) %>%
cols_align(
align = "left",
columns = statistic
) %>%
cols_align(
align = "right",
columns = attachment
) %>%
tab_options(row_group.font.weight = "bold",
row_group.background.color = brewer.pal(9,"Greys")[3],
table.font.size = 12,
data_row.padding = 3,
row_group.padding = 4,
summary_row.padding = 2,
column_labels.font.size = 14,
row_group.font.size = 12,
table.width = pct(100)) %>%
tab_header(
title = md("**summary of attachment weights**"), # Add a bold title
subtitle = md("*expressed as relative percentage of scaled mass index or raw body mass*") # Add an italic subtitle
)
```
#### sample size summary
```{r, echo=FALSE}
body_traits_smi %>%
group_by(tag_type) %>%
summarise(n = n_distinct(ring))
```
### distribution of time differences between repeated measures
```{r}
hist (as.numeric (body_traits_smi$ date_diff))
body_traits_smi %>% filter (! is.na (date_diff)) %>% pull (date_diff) %>% min ()
body_traits_smi %>% filter (! is.na (date_diff)) %>% pull (date_diff) %>% max ()
```
### modelling
```{r, eval=FALSE}
smi_wing_lmer_DOY_season <- lme4::lmer(
smi_wing ~ tag_type * measure + scale(day_of_year) + (1 | ring),
data = body_traits_smi, REML = TRUE
)
```
### contrast interactive and non-interactive model
```{r}
smi_wing_lmer_DOY_season_interaction <-
lmer (smi_wing ~ tag_type * measure + scale (day_of_year) + (1 | ring),
data = body_traits_smi, REML = FALSE )
smi_wing_lmer_DOY_season_noninteraction <-
lmer (smi_wing ~ tag_type + measure + scale (day_of_year) + (1 | ring),
data = body_traits_smi, REML = FALSE )
```
```{r, echo=FALSE}
anova_results <- anova(smi_wing_lmer_DOY_season_interaction,
smi_wing_lmer_DOY_season_noninteraction, test = "Chisq") %>%
broom::tidy()
anova_results %>%
kable("html", digits = 3, caption = "Likelihood Ratio Test for Model Comparison") %>%
kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))
```
### statistics
```{r, eval=FALSE}
# set seed for reproducibility
set.seed(1234)
# Derive confidence intervals of effect sizes from parametric bootstrapping
tidy_smi_wing_lmer_DOY_season <-
tidy(smi_wing_lmer_DOY_season, conf.int = TRUE, conf.method = "boot", nsim = 1000)
# run rptR to obtain repeatabilities of random effects
rpt_smi_wing_lmer_DOY_season <-
rpt(smi_wing ~ tag_type * measure + scale(day_of_year) + (1 | ring),
grname = c("ring", "Fixed"),
data = body_traits_smi,
datatype = "Gaussian",
nboot = 1000, npermut = 1000, ratio = TRUE,
adjusted = TRUE, ncores = 4, parallel = TRUE)
# run partR2 on each model to obtain marginal R2, parameter estimates, and beta
# weights
R2m_smi_wing_lmer_DOY_season <-
partR2(smi_wing_lmer_DOY_season,
partvars = c("tag_type",
"measure",
"scale(day_of_year)"),
R2_type = "marginal",
nboot = 1000,
CI = 0.95,
max_level = 1)
R2c_smi_wing_lmer_DOY_season <-
partR2(smi_wing_lmer_DOY_season,
partvars = c("tag_type",
"measure",
"scale(day_of_year)"),
R2_type = "conditional",
nboot = 1000,
CI = 0.95,
max_level = 1)
stats_smi_wing_lmer_DOY_season <-
list(mod = smi_wing_lmer_DOY_season,
tidy = tidy_smi_wing_lmer_DOY_season,
rptR = rpt_smi_wing_lmer_DOY_season,
partR2m = R2m_smi_wing_lmer_DOY_season,
partR2c = R2c_smi_wing_lmer_DOY_season,
data = body_traits_smi)
```
```{r, eval=FALSE}
# saveRDS(stats_smi_wing_lmer_DOY_season, here("out/smi_date_results.rds"))
```
```{r, echo = FALSE}
stats_smi_wing_lmer_DOY_season <- readRDS(here("out/smi_date_results.rds"))
#### Table of effect sizes ----
# Retrieve sample sizes
sample_sizes <-
stats_smi_wing_lmer_DOY_season$data %>%
ungroup() %>%
summarise(Individual = n_distinct(ring),
Observations = nrow(.))
sample_sizes <-
as.data.frame(t(as.data.frame(sample_sizes))) %>%
rownames_to_column("term") %>%
rename(estimate = V1) %>%
mutate(stat = "n")
# clean model component names
mod_comp_names <-
data.frame(comp_name = c("Intercept (control group 1st capture)",
"Julian date of measure",
"GPS group 1st capture",
"PTT group 1st capture",
"Control group 2nd capture (within season)",
"Control group 2nd capture (between seasons)",
"GPS group*2nd capture (within season)",
"PTT group*2nd capture (within season)",
"GPS group*2nd capture (between seasons)",
"PTT group*2nd capture (between seasons)",
"Total Marginal \U1D479\U00B2",
"Tracking device type",
"Measurement sequence",
"Julian date of measure",
"Total Conditional \U1D479\U00B2",
"Individual",
"Residual",
"Individual",
"Residual",
"Individuals",
"Observations"))
# Fixed effect sizes (non-standardized)
fixefTable <-
stats_smi_wing_lmer_DOY_season$tidy %>%
dplyr::filter(effect == "fixed") %>%
dplyr::select(term, estimate, conf.low, conf.high) %>%
as.data.frame() %>%
mutate(stat = "fixed")
# Fixed effect sizes (standardized)
fixef_bw_Table <-
stats_smi_wing_lmer_DOY_season$partR2m$BW %>%
as.data.frame() %>%
mutate(stat = "fixed_bw") %>%
rename(conf.low = CI_lower,
conf.high = CI_upper)
# Semi-partial R2 estimates
R2Table <-
bind_rows(stats_smi_wing_lmer_DOY_season$partR2m$R2,
stats_smi_wing_lmer_DOY_season$partR2c$R2[1,]) %>%
dplyr::select(term, estimate, CI_lower, CI_upper) %>%
as.data.frame() %>%
mutate(stat = "partR2") %>%
rename(conf.low = CI_lower,
conf.high = CI_upper)
# Random effects variances
ranefTable <-
stats_smi_wing_lmer_DOY_season$tidy %>%
dplyr::filter(effect == "ran_pars") %>%
dplyr::select(group, estimate, conf.low, conf.high) %>%
as.data.frame() %>%
mutate(stat = "rand") %>%
rename(term = group) %>%
mutate(estimate = estimate^2,
conf.high = conf.high^2,
conf.low = conf.low^2)
# Adjusted repeatabilities
coefRptTable <-
stats_smi_wing_lmer_DOY_season$rptR$R_boot %>%
dplyr::select(-Fixed) %>%
mutate(residual = 1 - rowSums(.)) %>%
apply(., 2,
function(x) c(mean (x), quantile (x, prob = c(0.025, 0.975)))) %>%
t() %>%
as.data.frame() %>%
rownames_to_column("term") %>%
rename(estimate = V1,
conf.low = `2.5%`,
conf.high = `97.5%`) %>%
mutate(stat = "RptR")
# Store all parameters into a single table and clean it up
allCoefs_mod <-
bind_rows(fixefTable[1,],
fixef_bw_Table,
R2Table,
ranefTable,
coefRptTable,
sample_sizes) %>%
slice(c(1,6,2,3,5,4,9,10,7,8,11:21)) %>%
bind_cols(.,
mod_comp_names) %>%
mutate(coefString = ifelse(!is.na(conf.low),
paste0("[",
round(conf.low, 2), ", ",
round(conf.high, 2), "]"),
NA),
effect = c(rep("Fixed effects \U1D6FD (standardized; units = \u011D)", nrow(fixefTable[1,])),
rep("Fixed effects \U1D6FD (standardized; units = \u011D)", nrow(fixef_bw_Table)),
rep("Partitioned \U1D479\U00B2", nrow(R2Table)),
rep("Random effects \U1D70E\U00B2", nrow(ranefTable)),
rep("Adjusted repeatability \U1D45F", nrow(coefRptTable)),
rep("Sample sizes \U1D45B", nrow(sample_sizes)))) %>%
dplyr::select(effect, everything())
# draw gt table
smi_wing_lmer_table <-
allCoefs_mod %>%
dplyr::select(effect, comp_name, estimate, coefString) %>%
gt(rowname_col = "row",
groupname_col = "effect") %>%
cols_label(comp_name = html("<i>Scaled mass index (SMI) dynamics</i>"),
estimate = "Mean estimate",
coefString = "95% confidence interval") %>%
fmt_number(columns = c(estimate),
rows = 1:19,
decimals = 2,
use_seps = FALSE) %>%
fmt_number(columns = c(estimate),
rows = 20:21,
decimals = 0,
use_seps = FALSE) %>%
sub_missing(columns = 1:4,
missing_text = "") %>%
cols_align(align = "left",
columns = c(comp_name)) %>%
tab_options(row_group.font.weight = "bold",
row_group.background.color = brewer.pal(9,"Greys")[3],
table.font.size = 12,
data_row.padding = 3,
row_group.padding = 4,
summary_row.padding = 2,
column_labels.font.size = 14,
row_group.font.size = 12,
table.width = pct(50))
smi_wing_lmer_table
```
```{r, eval=FALSE, include=FALSE, echo=FALSE}
gtsave(smi_wing_lmer_table %>% tab_options(table.width = pct(60)),
here("tabs_figs/table_1_.png"))
```
```{r}
gtsave (smi_wing_lmer_table, here ("tabs_figs/table_1.rtf" ))
```
# Breeding status and outcome
```{r}
breed_fates <-
read.csv (here ("data/breed.csv" ))
```
### Breeding data summary
```{r, echo = FALSE}
breed_fates %>%
dplyr::summarise(n_nests = n_distinct(nest_id)) %>%
kable(format = "html", caption = "number of distinct nests") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = FALSE)
breed_fates %>%
group_by(parent, tag_type) %>%
summarise(n_seasons = n_distinct(season)) %>%
filter(tag_type != "control") %>%
table()
breed_fates %>% filter(parent == "RBBR")
breed_fates %>%
dplyr::select(nest_id, state) %>%
distinct() %>%
group_by(state) %>%
summarise(n_nests = n()/117) %>%
kable(format = "html", caption = "proportion of nests associated with each fate") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = FALSE)
breed_fates %>%
group_by(state) %>%
dplyr::summarise(n_nests = n_distinct(nest_id)) %>%
kable(format = "html", caption = "number of distinct nests for each fate") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = FALSE)
breed_fates %>%
group_by(tag_type) %>%
dplyr::summarise(n_nests = n_distinct(nest_id)) %>%
kable(format = "html", caption = "number of distinct nests for each treatment") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = FALSE)
breed_fates %>%
group_by(tag_type) %>%
dplyr::summarise(n_ind = n_distinct(parent)) %>%
kable(format = "html", caption = "number of distinct birds for each treatment") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = FALSE)
breed_fates %>%
group_by(tag_type) %>%
dplyr::summarise(n_ind = n_distinct(parent)) %>%
kable(format = "html", caption = "number of distinct birds for each treatment") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = FALSE)
```
## check if breeding fates of subsequent seasons with the treatment are different than the deployment season
```{r, eval=FALSE}
# Determine first treatment season per parent and tag_type
first_treatment_season <- breed_fates %>%
filter(tag_type %in% c("GPS", "PTT")) %>%
group_by(parent, tag_type) %>%
summarise(first_treatment = min(season), .groups = "drop")
# Join back to main data and classify treatment stage
breed_fates_ <- breed_fates %>%
left_join(first_treatment_season, by = c("parent", "tag_type")) %>%
mutate(
treatment_stage = case_when(
tag_type == "control" ~ "control",
season == first_treatment ~ "first_treatment_season",
season > first_treatment ~ "subsequent_treatment_season",
TRUE ~ NA_character_
)
) %>%
mutate(treatment_stage = relevel(factor(treatment_stage),
ref = "control"))
# set seed for reproducibility
set.seed(1234)
breed_brm2 <- brm(
state ~ treatment_stage + (1 | season) + (1 | parent),
data = breed_fates_,
family = categorical(link = "logit"),
chains = 4,
iter = 4000,
control = list(adapt_delta = 0.99, max_treedepth = 15)
)
# consolidate results
stats_breed_brm_ <-
list(mod = breed_brm2,
data = breed_fates_)
```
```{r, echo=FALSE, results='hide'}
# saveRDS(stats_breed_brm_, here("out/breed_carryover_season.rds"))
stats_breed_brm_ <-
readRDS(here("out/breed_carryover_season.rds"))
predictions_treatment_stage <-
ggpredict(stats_breed_brm_$mod,
terms = "treatment_stage",
type = "fixed",
bias_correction = TRUE
) %>%
as.data.frame()
facet_labels <- c(
"aband" = "abandoned",
"brood" = "hatched",
"fail" = "failed",
"fledge" = "fledged"
)
# Assuming 'nest_fates' is your data frame
# Calculate the sample sizes
sample_sizes <- stats_breed_brm_$data %>%
group_by(treatment_stage, state) %>%
summarise(n = n(),
n_nests = n_distinct(nest_id)) %>%
ungroup() %>%
rename(x = treatment_stage,
response.level = state) %>%
arrange(x, response.level)
# Merge the sample sizes with your predictions
predictions_treatment_stage <-
predictions_treatment_stage %>%
left_join(sample_sizes, by = c("x", "response.level")) %>%
mutate(response.level = factor(response.level, levels = c("aband", "fail", "brood", "fledge")))
custom_colors_ <- brewer.pal(9, "Set1")[c(9, 3, 2)]
bdot_treatment_stage_facets <- c('control' = "control",
'first_treatment_season' = "deployment season",
'subsequent_treatment_season' = "subsequent season")
# treatment_stage_breeding_fate_plot_ <-
ggplot(data = predictions_treatment_stage) +
geom_errorbar(aes(x = response.level, y = predicted, ymin = conf.low, ymax = conf.high, color = x),
size = 0.3, linetype = "solid", width = 0.1) +
geom_point(aes(x = response.level, y = predicted, fill = x, shape = response.level), size = 4, color = "black", alpha = 1) +
# scale_fill_manual(values = custom_colors) +
# scale_color_manual(values = custom_colors_) +
geom_text(aes(label = n_nests, y = -0.2, x = response.level),
size = 3, color = "black", vjust = 0) +
theme(legend.position = "none",
text = element_text(family = "Verdana"),
axis.ticks = element_blank(),
axis.title.y = element_text(size = 11, colour = "black"),
axis.title.x = element_blank(),
axis.text.y = element_text(size = 10),
axis.text.x = element_text(size = 10, colour = "black", angle = 45, hjust = 1),
plot.margin = unit(c(0.2,0.2,0.2,0.2), "cm"),
panel.background = element_rect(fill = 'transparent', color = NA),
plot.background = element_rect(fill = 'transparent', color = NA),
strip.background = element_rect(fill = 'grey90', color = NA),
# strip.text = element_blank(),
panel.spacing = unit(1, "cm"),
panel.grid.major.y = element_line(size = 0.25, lineend = "round", colour = "grey90"),
panel.grid.minor = element_blank(),
panel.grid.major.x = element_blank()) +
ylab("fate probability\n± 95% CrI") +
xlab("breeding attempt fate") +
scale_y_continuous(limits = c(-0.2, 1), breaks = c(0, 0.25, 0.5, 0.75, 1)) +
scale_x_discrete(labels = c("abandoned", "failed", "hatched", "fledged")) +
facet_grid(~ x, labeller = as_labeller(bdot_treatment_stage_facets)) +
scale_shape_manual(values = c(21, 22, 23, 24, 25))
```
### modelling
```{r, eval=FALSE}
# set seed for reproducibility
set.seed(1234)
# Fit the mixed-effects multinomial logistic regression model
breed_brm <- brm(
state ~ tag_type + (1 | season) + (1 | parent),
data = breed_fates,
family = categorical(link = "logit"),
chains = 4,
iter = 4000,
control = list(adapt_delta = 0.99, max_treedepth = 15)
)
# consolidate results
stats_breed_brm <-
list(mod = breed_brm,
data = breed_fates)
```
# Apparent survival
### RMark data preparation
```{r, eval=FALSE}
ch <- read.csv(here("data/ch.csv"))
# Specify the states and create the dataframe for RMark
ch_processed <- process.data(ch, model = "Multistrata", groups = NULL)
ch_ddl <- make.design.data(ch_processed)
# Define parameter formulas
S.model <- list(formula = ~ stratum) # Survival depends on the state
p.model <- list(formula = ~ stratum) # Detection probability depends on the state
Psi.model.stratum <- list(formula = ~ -1 + stratum)
Psi.model.stratumto <- list(formula = ~ -1 + stratum:tostratum)
Psi.model.stratumtotime <- list(formula = ~ -1 + stratum : tostratum : time) # Transition probabilities depend on the state and destination state
```
### RMark modelling
```{r, eval=FALSE}
model.list <- create.model.list("Multistrata")
surv.results <- mark.wrapper(model.list, data = ch_processed, ddl = ch_ddl,
threads = 4, brief = TRUE, delete = TRUE)
```
```{r, eval=FALSE, echo=FALSE}
# extract and format the transformed parameter estimates
surv_estimates <-
surv.results$S.model.p.model.Psi.model.stratumtotime$results$real %>%
bind_cols(data.frame(str_split_fixed(rownames(.), " ",
n = 5)), .) %>%
dplyr::mutate(tag_type = as.factor(ifelse(unlist(str_extract_all(X2,"[CGP]")) == "C", "control",
ifelse(unlist(str_extract_all(X2,"[CGP]")) == "G", "GPS", "PTT"))),
parameter = as.factor(ifelse(X1 == "S", "Phi",
ifelse(X1 == "p", "p", "Psi"))),
transition = ifelse(X1 == "Psi", paste0(str_sub(X2, 2, 2), X3), NA),
year = ifelse(X1 == "Psi", ifelse(str_sub(X5, 10, 11) == "t1", 2021, 2022), NA)) %>%
dplyr::select(parameter, tag_type, transition, year, estimate, lcl, ucl, se) %>%
`rownames<-`( NULL )
```
# Figure 2 - multi-panelled plot of all results
```{r, echo=FALSE}
## scaled mass index plot
stats_smi_wing_lmer <- readRDS(here("out/smi_date_results.rds"))
sample_sizes <-
stats_smi_wing_lmer$data %>%
group_by(tag_type, measure) %>%
summarise(n = n()) %>%
ungroup() %>%
mutate(n = as.character(n)) %>%
mutate(n = ifelse(tag_type %in% c("GPS", "PTT") & measure == "second", paste0(n, "*"), n))
# repeated_measures_ <-
# stats_smi_wing_lmer$data %>%
# mutate(x_line = ifelse(measure == "first", as.numeric(factor(measure)) + 0.15,
# ifelse(measure == "second", as.numeric(factor(measure)) - 0.15, NA)),
# x_point = ifelse(measure == "first", as.numeric(factor(measure)),
# ifelse(measure == "second", as.numeric(factor(measure)), NA))) %>%
# left_join(sample_sizes, by = c("tag_type", "measure"))
repeated_measures_ <-
stats_smi_wing_lmer$data %>%
mutate(x_line = ifelse(measure == "first", 1 + 0.15,
ifelse(measure == "second_within",
2 - 0.15,
ifelse(measure == "second_between",
3 - 0.15, NA))),
x_point = ifelse(measure == "first", 1,
ifelse(measure == "second_within",
2,
ifelse(measure == "second_between",
3, NA)))) %>%
left_join(sample_sizes, by = c("tag_type", "measure"))
bdot_treatment_group_facets <- c('GPS' = "1.2g GPS tag",
'PTT' = "1.8 & 2.0g PTT tags",
'control' = "Control")
custom_colors <- brewer.pal(9, "Pastel1")[c(9, 3, 2)]
repeated_measures_plot <-
ggplot(data = repeated_measures_ %>% mutate(tag_type_ = ifelse(ring %in% c("CP16918", "CP9068", "CP16905"), "PTT2", tag_type))) +
geom_line(aes(x = x_line, y = smi_wing, group = ring, linetype = tag_type_), color = "black", alpha = 0.5) +
geom_jitter(size = 4, aes(x = x_point, y = smi_wing, fill = tag_type, shape = measure),
shape = 21, color = "black", alpha = 1, width = 0.075) +
scale_fill_manual(values = custom_colors) +
scale_linetype_manual(values = c("solid", "solid", "solid", "dashed")) +
facet_grid(~ tag_type, labeller = as_labeller(bdot_treatment_group_facets)) +
theme(legend.position = "none",
text = element_text(family = "Verdana"),
axis.title.x = element_blank(),
axis.ticks = element_blank(),
axis.text = element_text(size = 10, colour = "black"),
axis.title.y = element_text(size = 11, colour = "black"),
axis.text.y = element_text(size = 10, colour = "black"),
plot.margin = unit(c(0.2,0.2,0.2,0.2), "cm"),
panel.spacing = unit(1, "cm"),
panel.background = element_rect(fill = 'transparent', color = NA),
plot.background = element_rect(fill = 'transparent', color = NA),
strip.background = element_rect(fill = 'grey90', color = NA),
strip.text = element_text(size = 13, colour = "black"),
panel.grid.major.y = element_line(size = 0.25, lineend = "round", colour = "grey90"),
panel.grid.minor = element_blank(),
panel.grid.major.x = element_blank()) +
scale_x_continuous(limits = c(0.8, 3.2), breaks = c(1, 2, 3), labels = c("1st\ncapture", "2nd\nwithin\nseason", "2nd\nbetween\nseason")) +
scale_y_continuous(breaks = c(40, 50, 60, 70, 80), limits = c(37, 82)) +
xlab("capture event") +
ylab(expression("Scaled body mass index (" * hat(g) * ")")) +
geom_text(data = repeated_measures_ %>% ungroup() %>% dplyr::select(x_point, n, tag_type) %>% distinct(), aes(label = n, y = 37, x = x_point),
size = 3, color = "black", vjust = 0)
```
```{r, echo=FALSE, results='hide'}
## breeding fate plot
stats_breed_brm <-
readRDS(here("out/breed_results.rds"))
facet_labels <- c(
"aband" = "abandoned",
"brood" = "hatched",
"fail" = "failed",
"fledge" = "fledged"
)
# Generate predictions
predictions_tag_type_ <-
ggpredict(stats_breed_brm$mod,
terms = "tag_type",
type = "fixed",
bias_correction = TRUE
) %>%
as.data.frame()
# Assuming 'nest_fates' is your data frame
# Calculate the sample sizes
sample_sizes <- stats_breed_brm$data %>%
group_by(tag_type, state) %>%
summarise(n = n(),
n_nests = n_distinct(nest_id)) %>%
ungroup() %>%
rename(x = tag_type,
response.level = state) %>%
bind_rows(., data.frame(x = c("GPS", "GPS"), response.level = c("fledge", "aband"), n = c(0, 0), n_nests = c(0, 0))) %>%
arrange(x, response.level)
# Merge the sample sizes with your predictions
predictions_tag_type_ <-
predictions_tag_type_ %>%
left_join(sample_sizes, by = c("x", "response.level")) %>%
mutate(response.level = factor(response.level, levels = c("aband", "fail", "brood", "fledge")))
custom_colors_ <- brewer.pal(9, "Set1")[c(9, 3, 2)]
tag_type_breeding_fate_plot_ <-
ggplot(data = predictions_tag_type_) +
geom_errorbar(aes(x = response.level, y = predicted, ymin = conf.low, ymax = conf.high, color = x),
size = 0.3, linetype = "solid", width = 0.1) +
geom_point(aes(x = response.level, y = predicted, fill = x, shape = response.level), size = 4, color = "black", alpha = 1) +
scale_fill_manual(values = custom_colors) +
scale_color_manual(values = custom_colors_) +
geom_text(aes(label = n_nests, y = -0.2, x = response.level),
size = 3, color = "black", vjust = 0) +
theme(legend.position = "none",
text = element_text(family = "Verdana"),
axis.ticks = element_blank(),
axis.title.y = element_text(size = 11, colour = "black"),
axis.title.x = element_blank(),
axis.text.y = element_text(size = 10),
axis.text.x = element_text(size = 10, colour = "black", angle = 45, hjust = 1),
plot.margin = unit(c(0.2,0.2,0.2,0.2), "cm"),
panel.background = element_rect(fill = 'transparent', color = NA),
plot.background = element_rect(fill = 'transparent', color = NA),
strip.background = element_rect(fill = 'transparent', color = NA),
strip.text = element_blank(),
panel.spacing = unit(1, "cm"),
panel.grid.major.y = element_line(size = 0.25, lineend = "round", colour = "grey90"),
panel.grid.minor = element_blank(),
panel.grid.major.x = element_blank()) +
ylab("Fate probability\n± 95% CrI") +
xlab("Breeding attempt fate") +
scale_y_continuous(limits = c(-0.2, 1), breaks = c(0, 0.25, 0.5, 0.75, 1)) +
scale_x_discrete(labels = c("Abandoned", "Failed", "Hatched", "Fledged")) +
facet_grid(~ x, labeller = as_labeller(bdot_treatment_group_facets)) +
scale_shape_manual(values = c(21, 22, 23, 24, 25))
```
```{r, echo = FALSE}
## apparent survival plot
surv_estimates <- readRDS(here("out/survival_results.rds"))
ch <- read.csv(here("data/ch.csv"))
sample_sizes_surv <-
ch %>%
mutate(tag_type = ifelse(str_detect(ch, "P"), "PTT",
ifelse(str_detect(ch, "G"), "GPS",
ifelse(str_detect(ch, "C"), "control", NA)))) %>%
group_by(tag_type) %>%
summarise(n = n())
tag_type_surv_p_plot <-
ggplot(data = filter(surv_estimates, parameter %in% c("Phi", "p")),
aes(x = parameter, y = estimate, fill = tag_type)) +
geom_errorbar(aes(ymin = lcl, ymax = ucl, color = tag_type),
size = 0.3, linetype = "solid", width = 0.05) +
geom_point(aes(shape = parameter), size = 4, color = "black", alpha = 1) +
scale_fill_manual(values = custom_colors) +
scale_color_manual(values = custom_colors_) +
theme(legend.position = "none",
text = element_text(family = "Verdana"),
axis.title.x = element_blank(),
axis.ticks = element_blank(),
axis.title.y = element_text(size = 11, colour = "black"),
axis.text.y = element_text(size = 10),
axis.text.x = element_text(size = 10, colour = "black"),
plot.margin = unit(c(0.2,0.2,0.2,0.2), "cm"),
panel.spacing = unit(1, "cm"),
panel.background = element_rect(fill = 'white', color = NA),
plot.background = element_rect(fill = 'white', color = NA),
strip.background = element_rect(fill = 'white', color = NA),
strip.text = element_blank(),
panel.grid.major.y = element_line(size = 0.25, lineend = "round", colour = "grey90"),
panel.grid.minor = element_blank(),
panel.grid.major.x = element_blank()) +
ylab("Probability\n± 95% CI") +
scale_y_continuous(limits = c(-0.15, 1.1), breaks = c(0, .25, .5, .75, 1)) +
scale_x_discrete(labels = c("\u03C6", expression(italic(p)))) +
scale_shape_manual(values = c(21, 22)) +
facet_grid(~ tag_type, labeller = as_labeller(bdot_treatment_group_facets)) +
geom_text(data = sample_sizes_surv, aes(x = 1.5, y = -0.15, label = n), size = 3, color = "black", vjust = 0)
```
```{r, echo=FALSE, fig.height=10, fig.width=8}
combo_plot <-
repeated_measures_plot / plot_spacer() / tag_type_breeding_fate_plot_ / tag_type_surv_p_plot +
plot_layout(heights = c(2, 0.1, 1, 1)) +
plot_annotation(
tag_levels = 'a',
tag_prefix = "(",
tag_suffix = ")",
theme = theme(plot.tag = element_text(family = "Verdana", size = 14, face = "bold"))
)
combo_plot
ggsave(plot = combo_plot,
filename = "Figure_2.svg",
path = here("tabs_figs/"), dpi = 300,
height = 10, width = 8, units = "in")
```
# Behavioural Observations
### data summary
```{r, echo=FALSE}
behav <-
read.csv(here("data/behav.csv")) %>%
mutate(time = as_hms(time))
ggplot(behav %>% filter(!is.na(time)), aes(x = time)) +
geom_histogram(binwidth = 3600, fill = "skyblue", color = "black") +
scale_x_time(breaks = as_hms(seq(0, 86399, by = 14400)),
labels = time_format("%H:%M")) +
labs(title = "diurnal variation in the time of behavioural observations", x = "Time of Day", y = "Count") +
theme_minimal()
behav %>%
group_by(tagged) %>%
summarise(n_obs = n()) %>%
mutate(tagged = ifelse(tagged == "Y", "tagged", "untagged")) %>%
kable(format = "html", caption = "summary of behavioural observations") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = FALSE)
```
### energy budget estimation
```{r, echo=FALSE}
# Table 5 from:
# Morrier, A. and McNeil, R., 1991. Time-activity budget of Wilson's and Semipalmated Plovers in a tropical environment. The Wilson Bulletin, pp.598-620.
Morrier_McNeil_table5 <-
data.frame(
Sampling_Sessions = c("18-21 October", "11-12 November", "21-22 November", "1 December",
"14-16 January", "29-30 January", "15-17 February", "8-9 March",
"22-24 March", "11-13 April"),
Feeding = c(18.38, 3.49, 5.01, 3.28, 2.18, 2.30, 0.39, 3.72, 3.70, 14.83),
Preening = c(0.58, 0.29, 0.73, 0.35, 0.86, 0.77, 0.54, 0.24, 0.64, 0.39),
Resting = c(8.56, 17.44, 15.89, 17.53, 17.58, 17.53, 19.60, 17.70, 17.81, 10.83),
Flight = c(0.34, 1.47, 0.67, 0.02, 0.56, 0.32, 0.00, 0.79, 0.15, 1.02),
Walking = c(0.12, 0.18, 0.62, 0.31, 0.35, 0.57, 0.09, 0.06, 0.17, 0.38),
Fight = c(0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.06, 0.21),
Ground_Display = c(0.00, 0.07, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.58, 0.03),
Alert = c(0.62, 0.75, 1.03, 0.48, 0.80, 1.05, 0.75, 1.46, 0.75, 1.73),
Running = c(0.03, 0.12, 0.03, 0.00, 0.11, 0.09, 0.00, 0.12, 0.04, 0.32),
Nocturnal_Resting = c(21.18, 21.39, 21.42, 21.72, 21.65, 21.50, 21.26, 20.96, 20.70, 20.43),
Total_Energy_Expenditure = c(49.81, 45.20, 45.43, 43.70, 44.08, 44.10, 42.63, 45.06, 44.61, 50.17)
)
Morrier_McNeil_table5 %>%
mutate(across(
.cols = c(Feeding, Preening, Resting, Flight, Walking, Fight,
Ground_Display, Alert, Running, Nocturnal_Resting),
.fns = ~ (. / Total_Energy_Expenditure) * 100,
.names = "prop_{.col}"
)) %>%
summarise(across(
.cols = starts_with("prop_"),
.fns = mean,
.names = "avg_{.col}"
)) %>% t() %>%
kable(format = "html", caption = "Table 5: Morrier and McNeil (1991). Time-activity budgets of Wilson's and Semipalmated Plovers in a tropical environment") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = FALSE)
```
### modelling
```{r}
inc_brood_model <-
glmer (inc_brood ~ tagged + sin_time + cos_time + (1 | bird),
data = behav, family = binomial)
```
```{r, echo=FALSE}
# Get tidy output
tidy_table_i <- tidy(inc_brood_model, effects = "fixed")
# Render as a formatted table
tidy_table_i %>%
kable("html", digits = 3, caption = "Fixed Effects of the GLMM on incubation & brooding behaviour") %>%
kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))
```
```{r}
walk_model <-
glmer (walk ~ tagged + sin_time + cos_time + (1 | bird),
data = behav, family = binomial)
```
```{r, echo=FALSE}
# Get tidy output
tidy_table_w <- tidy(walk_model, effects = "fixed")
# Render as a formatted table
tidy_table_w %>%
kable("html", digits = 3, caption = "Fixed Effects of the GLMM on walking behaviour") %>%
kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))
```
```{r}
feed_model <-
glmer (feed ~ tagged + sin_time + cos_time + (1 | bird),
data = behav, family = binomial)
```
```{r, echo=FALSE}
# Get tidy output
tidy_table_f <- tidy(feed_model, effects = "fixed")
# Render as a formatted table
tidy_table_f %>%
kable("html", digits = 3, caption = "Fixed Effects of the GLMM on feeding behaviour") %>%
kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))
```
```{r}
fly_model <-
glmer (fly ~ tagged + sin_time + cos_time + (1 | bird),
data = behav, family = binomial)
```
```{r, echo=FALSE}
# Get tidy output
tidy_table_fl <- tidy(fly_model, effects = "fixed")
# Render as a formatted table
tidy_table_fl %>%
kable("html", digits = 3, caption = "Fixed Effects of the GLMM on flying behaviour") %>%
kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))
```
# Movements and tag operation
```{r, echo=FALSE}
tag_metadata <- read.csv(here("data/tag_metadata.csv"))
```
```{r, echo=FALSE}
tag_metadata %>%
filter(tag_type != "control") %>%
group_by(tag_type) %>%
summarise(mean = mean(deploy_duration, na.rm = TRUE)/7,
max = max(deploy_duration, na.rm = TRUE)/7,
min = min(deploy_duration, na.rm = TRUE)/7) %>%
kable(format = "html", caption = "summary of deployment duration (weeks)") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = FALSE)
```
```{r, echo=FALSE}
move <- read.csv(here("data/move.csv"))
```
```{r, echo=FALSE}
move %>%
group_by(parent, tag_type) %>%
summarise(total_fixes = n(),
max_time = max(as.numeric(time_since_deployment), na.rm = TRUE)) %>%
rowwise() %>%
mutate(max_time_weeks = (max_time / 60 / 60 / 24 / 7)) %>%
arrange(max_time) %>%
mutate(fixes_per_week = total_fixes/max_time_weeks) %>%
group_by(tag_type) %>%
summarise(mean = mean(fixes_per_week),
sd = sd(fixes_per_week),
min = min(fixes_per_week),
max = max(fixes_per_week)) %>% t() %>%
kable(format = "html", caption = "summary of fixes obtained per week") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = FALSE)
```
```{r, echo=FALSE}
move %>%
group_by(tag_type) %>%
summarise(total_points = n()) %>%
kable(format = "html", caption = "summary of total number of fixes acquired") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = FALSE)
```
```{r, echo=FALSE}
move %>%
group_by(parent, tag_type) %>%
summarise(total_fixes = n(),
max_time = max(time_since_deployment, na.rm = TRUE)) %>%
group_by(tag_type) %>%
summarise(mean = mean(max_time)/86400/7,
sd = sd(max_time)/86400/7,
min = min(max_time)/86400/7,
max = max(max_time)/86400/7) %>%
kable(format = "html", caption = "summary of operational duration (weeks)") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = FALSE)
```
### plot of tag data collected over time
```{r, echo=FALSE}
custom_fill_ <- brewer.pal(9, "Pastel1")[c(2, 3)]
tag_type_facets <- c('GPS' = "1.2g archival GPS tag",
'PTT' = "1.8g & 2.0g* Argos PTT tags")
facets <-
data.frame(parent = c("RL-YL", "RW-BG", "RL-YR", "RR-BR", "RB-BY", "RW-BY", "RB-BR", "RR-LR", "RR-YY", "RW-LB", rev(c("RB-OG", "RB-GB", "RR-OW", "RB-WO", "RR-GO", "RR-RG", "RB-OY", "RB-WR", "RB-RL", "RW-LG"))),
label = c("tag retrieved but no data onboard", "tag retrieved but no data onboard", "bird did not return", "bird did not return", rep(NA, 16)),
tag_type = c(rep("GPS", 10), rep("PTT", 10))) %>%
mutate(parent = factor(parent,
levels = c("RL-YL", "RW-BG", "RL-YR", "RR-BR", "RB-BY", "RW-BY", "RB-BR", "RR-LR", "RR-YY", "RW-LB", rev(c("RB-OG", "RB-GB", "RR-OW", "RB-WO", "RR-GO", "RR-RG", "RB-OY", "RB-WR", "RB-RL", "RW-LG")))))
# Define the bin size
bin_width_ <- 30/7
# Bin the data
binned_data <-
move %>%
mutate(ring = factor(ring, levels = c("CP9098" , "CP9042" , "CP9055" , "CP9032" , "CP9096" , "CP9066" , "CP16134", "CP16111", "CP16607" ,"CP16617", "CP16918", "CP16905" ,"CP16609", "CP16616", "CP9068" , "CP9099"))) %>%
mutate(parent = factor(parent,
levels = c("RLYL", "RWBG", "RLYR", "RRBR", "RBBY", "RWBY", "RBBR", "RRLR", "RRYY", "RWLB", rev(c("RBOG", "RBGB", "RROW", "RBWO", "RRGO", "RRRG", "RBOY", "RBWR", "RBRL", "RWLG"))))) %>%
mutate(time_days = as.numeric(time_since_deployment/86400),
time_weeks = as.numeric(time_since_deployment/604800)) %>%
mutate(
x_bin_ = cut(time_weeks, breaks = seq(floor(min(time_weeks)), ceiling(max(time_weeks)), by = bin_width_))) %>%
group_by(parent, x_bin_) %>%
summarise(count = n(), .groups = 'drop',
month_ = round(mean(month(date_local)))) %>%
mutate(
x_bin_center = as.numeric(sub("\\((.+),.*", "\\1", x_bin_)) + bin_width_ / 2
) %>%
left_join(., tag_metadata %>% dplyr::select(code, tag_type) %>% rename(parent = code), by = "parent") %>%
distinct() %>%
group_by(x_bin_center) %>%
mutate(month__ = month(round(mean(month_)), label = TRUE)) %>%
ungroup()
custom_labels <- function(x) {
replacements <- c("RB-RL" = "*RB-RL",
"RR-RG" = "*RR-RG",
"RR-GO" = "*RR-GO")
ifelse(x %in% names(replacements), replacements[x], x)
}
tag_data_plot <-
binned_data %>%
mutate(parent = paste(str_sub(parent, 1, 2), str_sub(parent, 3, 4), sep = "-")) %>%
mutate(parent = factor(parent,
levels = rev(c("RL-YL", "RW-BG", "RL-YR", "RR-BR", "RB-BY", "RW-BY", "RB-BR", "RR-LR", "RR-YY", "RW-LB", rev(c("RB-OG", "RB-GB", "RR-OW", "RB-WO", "RR-GO", "RR-RG", "RB-OY", "RB-WR", "RB-RL", "RW-LG"))))),
tag_type = factor(as.factor(tag_type), levels = c("PTT", "GPS"))) %>%
# mutate(parent = factor(parent,
# levels = rev(c("RLYL", "RWBG", "RLYR", "RRBR", "RBBY", "RWBY", "RBBR", "RRLR", "RRYY", "RWLB", rev(c("RBOG", "RBGB", "RROW", "RBWO", "RRGO", "RRRG", "RBOY", "RBWR", "RBRL", "RWLG"))))),
# tag_type = factor(as.factor(tag_type), levels = c("PTT", "GPS"))) %>%
ggplot(aes(x = x_bin_center, y = parent)) +
geom_vline(xintercept = 0, color = "black", size = 0.2) +
geom_point(aes(size = count, fill = tag_type), shape = 21, alpha = 1) +
scale_size_continuous(range = c(1, 10),
breaks = c(1, 10, 100, 500),
labels = c("1", "10", "100", "500"), name = "fixes per month") +
theme_bw() +
theme(
text = element_text(family = "Verdana"),
legend.title = element_text(size = 10),
legend.text = element_text(size = 8),
strip.text = element_text(size = 10),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.ticks = element_line(linewidth = 0.3, colour = "grey40"),
axis.ticks.length = unit(0.2, "cm"),
panel.border = element_rect(linetype = "solid", colour = "grey"),
legend.position.inside = c(0.1, 0.9)
) +
geom_text(data = facets %>% mutate(tag_type = factor(as.factor(tag_type), levels = c("PTT", "GPS"))), aes(x = 1, y = parent, label = label), hjust = 0, size = 3, fontface = "italic") +
facet_wrap(ncol = 1, "tag_type", scales = "free_y", labeller = as_labeller(tag_type_facets), strip.position = "left") +
scale_fill_manual(values = custom_fill_) +
theme(legend.position = c(0.88, 0.14),
strip.text = element_text(size = 12),
strip.placement = "outside",
legend.box.background = element_rect(fill = "white", color = "grey50"),
text = element_text(family = "Verdana"),
strip.background = element_rect(fill = "white", colour = "white"),
panel.grid.major.y = element_line(size = 0.25, lineend = "round", colour = "grey90"),
axis.title = element_text(size = 11, colour = "black"),
axis.text = element_text(size = 10, colour = "black"),
axis.ticks.y = element_blank(),
axis.title.x.top = element_blank(),
axis.text.x.top = element_text(angle = 45, hjust = 0, vjust = 0)
) +
ylab("bird identity") +
xlab("weeks since tag deployment") +
guides(fill = "none", size = guide_legend(override.aes = list(color = "black", fill = "grey90"),
title.position = "top", title.hjust = 0.5,
label.position = "bottom", nrow = 1)) +
scale_x_continuous(expand = expansion(add = c(10, 3)),
sec.axis = sec_axis(~ ., name = "calander month",
breaks = binned_data %>% select(x_bin_center, month__) %>% distinct() %>%
slice(seq(1, n(), by = 4)) %>% pull(x_bin_center),
labels = binned_data %>% select(x_bin_center, month__) %>% distinct() %>%
slice(seq(1, n(), by = 4)) %>% pull(month__))) +
scale_y_discrete(labels = custom_labels)
```
```{r, echo=FALSE, warning=FALSE, message=FALSE, results='hide'}
move_sf <-
move %>%
filter(longitude > 100 & latitude < 0) %>%
filter(tag_type == "GPS" | (tag_type == "PTT" & locationClass %in% c(1, 2, 3))) %>%
filter(site == "KK") %>%
st_as_sf(.,
coords = c("longitude", "latitude"),
crs = "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
### Define the Lambert Conformal Conic projection parameters
lcc_params <- st_crs("+proj=lcc +lat_1=-42.41 +lat_2=-42.41 +lat_0=-42.41 +lon_0=173.68 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs")
### RBRL dashed track
RBRL_mokau_lon <-
move_sf %>%
filter(parent == "RBRL") %>%
filter(date_local >= as.POSIXct("2022-03-08 17:14:36", tz = "Pacific/Auckland")) %>%
filter(date_local <= as.POSIXct("2022-03-14 08:53:39", tz = "Pacific/Auckland")) %>%
sfc_as_cols(., names = c("loc_longitude", "loc_latitude")) %>%
st_drop_geometry() %>%
arrange(date_local) %>%
slice(1) %>% pull(loc_longitude)
RBRL_mokau_lat <-
move_sf %>%
filter(parent == "RBRL") %>%
filter(date_local >= as.POSIXct("2022-03-08 17:14:36", tz = "Pacific/Auckland")) %>%
filter(date_local <= as.POSIXct("2022-03-14 08:53:39", tz = "Pacific/Auckland")) %>%
sfc_as_cols(., names = c("loc_longitude", "loc_latitude")) %>%
st_drop_geometry() %>%
arrange(date_local) %>%
slice(1) %>% pull(loc_latitude)
RBRL_track <-
move_sf %>%
filter(parent == "RBRL") %>%
filter(date_local >= as.POSIXct("2022-03-07 17:14:36", tz = "Pacific/Auckland")) %>%
filter(date_local <= as.POSIXct("2022-03-15 08:53:39", tz = "Pacific/Auckland")) %>%
st_transform(lcc_params) %>%
sfc_as_cols(., names = c("loc_longitude", "loc_latitude")) %>%
arrange(date_local) %>%
dplyr::select(parent, tagID, loc_longitude, loc_latitude, date_local) %>%
distinct()
### Process movement data to determine maximum displacement per individual
max_dist_location <- move_sf %>%
sfc_as_cols(names = c("loc_longitude", "loc_latitude")) %>%
st_drop_geometry() %>%
filter(locationClass %in% c(NA, 1, 2, 3)) %>% # Retain high-quality location fixes
dplyr::select(parent, tagID, loc_longitude, loc_latitude, date_local) %>%
mutate(cap_longitude = 173.68, cap_latitude = -42.41) %>% # Define capture location
distinct() %>%
rowwise() %>%
mutate(dist_from_deploy = distHaversine(
p1 = matrix(c(loc_longitude, loc_latitude), ncol = 2),
p2 = matrix(c(cap_longitude, cap_latitude), ncol = 2)
) / 1000) %>% # Compute distance in km
group_by(parent) %>%
mutate(max_dist = max(dist_from_deploy)) %>%
filter(max_dist == dist_from_deploy) %>% # Retain max displacement per individual
distinct() %>%
arrange(desc(max_dist)) %>%
mutate(
loc_longitude = ifelse(parent == "RBRL", RBRL_mokau_lon, loc_longitude),
loc_latitude = ifelse(parent == "RBRL", RBRL_mokau_lat, loc_latitude)
)
### Prepare data for plotting movement trajectories
BADO_curve_plot_df <- max_dist_location %>%
filter(max_dist > 100) %>% # Retain movements > 100 km
dplyr::select(parent, date_local, cap_longitude, cap_latitude) %>%
rename(longitude = cap_longitude, latitude = cap_latitude) %>%
mutate(data_type = "capture") %>%
bind_rows(
max_dist_location %>%
dplyr::select(parent, date_local, loc_longitude, loc_latitude) %>%
rename(longitude = loc_longitude, latitude = loc_latitude) %>%
mutate(data_type = "movement")
) %>%
st_as_sf(coords = c("longitude", "latitude"), crs = "EPSG:4326")
### Define the reference point (Kaikoura)
Kaikoura <- data.frame(longitude = 173.68, latitude = -42.41) %>%
st_as_sf(coords = c("longitude", "latitude"), crs = "EPSG:4326") %>%
st_transform(lcc_params)
### Retrieve and crop political boundaries
nz_bound <- gisco_get_countries(country = "New Zealand", resolution = "03")
nz_bound_crop <- st_crop(nz_bound, c(xmin = 165, ymin = -48, xmax = 179, ymax = -33)) %>%
st_transform(lcc_params)
### Read and transform water bodies (rivers, lakes, etc.)
nz_water <-
st_read(here("data/NZL_water_areas_dcw.shp")) %>%
filter(HYC_DESCRI == "Perennial/Permanent") %>%
st_transform(lcc_params)
### Retrieve and project digital elevation model (DEM)
nz_dem <- elevatr::get_elev_raster(locations = nz_bound_crop, z = 6, clip = "locations")
nz_dem_proj <- terra::project(terra::rast(nz_dem), terra::crs(nz_bound_crop))
mdtdf <- as.data.frame(nz_dem_proj, xy = TRUE) %>% rename(alt = 3)
### Compute hillshade for terrain visualization
sl <- terra::terrain(nz_dem_proj, "slope", unit = "radians")
asp <- terra::terrain(nz_dem_proj, "aspect", unit = "radians")
hillmulti <- map(c(270, 15, 60, 330), ~terra::shade(sl, asp, angle = 45, direction = .x, normalize = TRUE)) %>%
terra::rast() %>% sum()
hillmultidf <- as.data.frame(hillmulti, xy = TRUE)
### Define movement trajectories
bbox <- st_bbox(BADO_curve_plot_df)
parent_groups <- BADO_curve_plot_df %>% group_by(parent) %>% mutate(n_points = n()) %>% filter(n_points > 1)
curved_lines <- list()
for (parent in unique(parent_groups$parent)) {
parent_data <- filter(parent_groups, parent == !!parent) %>% arrange(date_local)
for (i in 1:(nrow(parent_data) - 1)) {
start <- st_coordinates(parent_data$geometry[i])
end <- st_coordinates(parent_data$geometry[i + 1])
curved_lines[[length(curved_lines) + 1]] <- data.frame(x = c(start[1], end[1]), y = c(start[2], end[2]), parent = parent)
}
}
curved_lines_df <- do.call(rbind, curved_lines) %>%
st_as_sf(coords = c("x", "y"), crs = "EPSG:4326") %>%
st_transform(lcc_params) %>%
sfc_as_cols(names = c("t_x", "t_y")) %>%
st_drop_geometry()
### Extract movement curves for each individual
RWLG_curve <- curved_lines_df %>% filter(parent == "RWLG")
RBBR_curve <- curved_lines_df %>% filter(parent == "RBBR")
RRRG_curve <- curved_lines_df %>% filter(parent == "RRRG")
RBRL_curve <- curved_lines_df %>% filter(parent == "RBRL")
RBOY_curve <- curved_lines_df %>% filter(parent == "RBOY")
### Load and transform graticules for mapping
graticules <- st_graticule(ndiscr = 50) %>% st_transform(lcc_params) %>% st_geometry()
### Define color palette for elevation visualization
nz_pal <- hypsometric_tints_db %>%
filter(pal == "wiki-2.0_hypso") %>%
mutate(limit = seq(0, 3035, length.out = nrow(.))) %>%
mutate(hex = ifelse(limit > 1500, "#FFFFFF", hex))
```
```{r, echo=FALSE, fig.width=10, fig.height=10}
m <-
ggplot() +
geom_sf(data = nz_bound_crop, fill = "white", color = "white") +
list(
geom_raster(data = hillmultidf,
aes(x, y, fill = sum),
show.legend = FALSE,
alpha = .5),
scale_fill_distiller(palette = "Greys"),
new_scale_fill(),
geom_raster(data = mdtdf,
aes(x, y, fill = alt),
alpha = 0.5),
scale_fill_gradientn(
colors = nz_pal$hex,
values = scales::rescale(nz_pal$limit),
limit = range(nz_pal$limit)
)) %>%
blend("multiply") +
geom_sf(data = nz_water, color = NA, fill = "#698ecf", alpha = 0.85) +
geom_curve(data = RWLG_curve, aes(x = t_x[1], y = t_y[1], xend = t_x[2], yend = t_y[2]),
curvature = 0.2, linewidth = 0.5, color = "grey40") +
geom_curve(data = RBBR_curve, aes(x = t_x[1], y = t_y[1], xend = t_x[2], yend = t_y[2]),
curvature = 0.4, linewidth = 0.5, color = "grey40") +
geom_curve(data = RRRG_curve, aes(x = t_x[1], y = t_y[1], xend = t_x[2], yend = t_y[2]),
curvature = 0.40, linewidth = 0.5, color = "grey40") +
geom_curve(data = RBRL_curve, aes(x = t_x[1], y = t_y[1], xend = t_x[2], yend = t_y[2]),
curvature = 0.2, linewidth = 0.5, color = "grey40") +
geom_curve(data = RBOY_curve, aes(x = t_x[1], y = t_y[1], xend = t_x[2], yend = t_y[2]),
curvature = 0.2, linewidth = 0.5, color = "grey40") +
geom_point(data = BADO_curve_plot_df %>% filter(parent %in% c("RWLG", "RBBR", "RRRG", "RBOY")) %>% group_by(parent) %>% arrange(date_local) %>% slice(2),
aes(x = st_coordinates(geometry)[, 1], y = st_coordinates(geometry)[, 2], group = parent),
size = 3, shape = 21, fill = "grey40", stroke = 0.8, color = "grey40") +
geom_point(data = Kaikoura, aes(x = st_coordinates(geometry)[, 1], y = st_coordinates(geometry)[, 2]),
size = 5, shape = 21, fill = "grey40", stroke = 0.8, color = "grey40") +
geom_path(data = RBRL_track, aes(x = st_coordinates(geometry)[, 1], y = st_coordinates(geometry)[, 2]),
linewidth = 0.5, linetype = "dashed", color = "grey40") +
geom_point(data = RBRL_track %>% slice(n()), aes(x = st_coordinates(geometry)[, 1], y = st_coordinates(geometry)[, 2]),
size = 3, shape = 21, fill = "grey40", stroke = 0.8, color = "grey40") +
guides(fill = guide_colorsteps(barwidth = 20,
barheight = .5,
title.position = "right")) +
labs(fill = "m", x = "Longitude", y = "Latitude") +
coord_sf(xlim = c(bbox["xmin"] - 650000, bbox["xmax"] + 450000),
ylim = c(bbox["ymin"] - 550000, bbox["ymax"] + 850000),
crs = lcc_params) +
luke_theme +
theme(legend.position = "none",
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.ticks = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
panel.border = element_blank(),
panel.background = element_rect(fill = 'transparent', color = NA),
plot.background = element_rect(color = NA, fill = "transparent")) +
geom_text_repel(data = BADO_curve_plot_df %>%
filter(parent %in% c("RWLG", "RBBR", "RBOY")) %>%
group_by(parent) %>%
arrange(date_local) %>%
slice(2),
aes(x = st_coordinates(geometry)[, 1], y = st_coordinates(geometry)[, 2], group = parent, label = parent), force = 5, nudge_x = 100000, nudge_y = 10000) +
geom_text_repel(data = Kaikoura,
aes(x = st_coordinates(geometry)[, 1], y = st_coordinates(geometry)[, 2], group = parent), label = "Kaikōura", force = 5, nudge_x = 100000) +
geom_text_repel(data = BADO_curve_plot_df %>%
filter(parent %in% c("RRRG")) %>%
group_by(parent) %>%
arrange(date_local) %>%
slice(2),
aes(x = st_coordinates(geometry)[, 1], y = st_coordinates(geometry)[, 2], group = parent, label = parent), force = 5, nudge_x = -100000) +
geom_text_repel(data = RBRL_track %>% slice(n()),
aes(x = st_coordinates(geometry)[, 1], y = st_coordinates(geometry)[, 2], group = parent, label = parent), force = 5, nudge_x = -100000)
print(tag_data_plot + theme(plot.margin = margin(10, 100, 10, 10)))
# Overlay inset plot at specific location
vp <- viewport(x = 1.05, y = 0.85, width = 0.8, height = 0.85, just = c("right", "top"))
print(m, vp = vp)
```
# Supplementary Analysis
## Breeding status and outcome (pairs)
```{r, echo=FALSE, warning=FALSE, message=FALSE}
breed_fates_wide <-
breed_fates %>%
group_by(nest_id) %>%
mutate(row_number = row_number()) %>%
ungroup() %>%
pivot_wider(
id_cols = c(nest_id, season, state, date),
names_from = row_number,
values_from = c(parent, tag_type),
names_glue = "{.value}_{row_number}"
) %>%
mutate(
parent_1 = ifelse(is.na(parent_1), parent_2, parent_1),
tag_type_1 = ifelse(is.na(tag_type_1), as.character(tag_type_2), as.character(tag_type_1)),
parent_2 = ifelse(parent_1 == parent_2, NA, parent_2),
tag_type_2 = ifelse(parent_1 == parent_2, NA, as.character(tag_type_2))
) %>%
mutate(
swap = tag_type_1 == "control" & tag_type_2 %in% c("PTT", "GPS"),
parent_1 = ifelse(swap, parent_2, parent_1),
parent_2 = ifelse(swap, parent_1, parent_2),
tag_type_1 = ifelse(swap, tag_type_2, tag_type_1),
tag_type_2 = ifelse(swap, tag_type_1, tag_type_2)
) %>%
select(-swap) %>%
mutate(parent_2 = ifelse(is.na(parent_2), "unkn", parent_2),
tag_type_2 = ifelse(is.na(tag_type_2), "control", tag_type_2)) %>%
mutate(parent_combo = paste(parent_1, parent_2, sep = "-"),
tag_type_combo = paste(tag_type_1, tag_type_2, sep = "-"))
```
#### summary of pair breeding fate data
number of mates and tag-type treatments over the three-year period
```{r}
breed_fates_wide %>%
filter (parent_2 != "unkn" ) %>%
group_by (parent_1) %>%
summarise (n_mates = n_distinct (parent_2),
n_combos = n_distinct (tag_type_combo)) %>%
arrange (desc (n_mates))
```
number distinct nesting attempts with different treatment combinations
```{r}
breed_fates_wide %>%
group_by (tag_type_combo) %>%
summarise (n_nests = n_distinct (nest_id))
```
proportion of nests of different fates
```{r}
breed_fates_wide %>%
dplyr:: select (nest_id, state) %>%
distinct () %>%
group_by (state) %>%
summarise (n_nests = n ()/ 117 )
```
### modelling
```{r, eval=FALSE}
breed_pair_brm_ <- brm(
state ~ tag_type_combo + (1 | season),
data = breed_fates_wide,
family = categorical(link = "logit"),
chains = 4,
iter = 4000,
control = list(adapt_delta = 0.99, max_treedepth = 15)
)
stats_breed_pair_brm <-
list(mod = breed_pair_brm_,
data = breed_fates_wide)
```
```{r, echo=FALSE}
# saveRDS(stats_breed_pair_brm, here("out/breed_pair_results.rds"))
stats_breed_pair_brm <-
readRDS(here("out/breed_pair_results.rds"))
```
#### supplementary figure 1
```{r, echo=FALSE, warning=FALSE, message=FALSE, error=FALSE, cache=FALSE, results='hide'}
# Generate predictions
predictions_tag_type_pair <-
suppressMessages(suppressWarnings(
ggpredict(stats_breed_pair_brm$mod,
terms = "tag_type_combo",
type = "fixed",
bias_correction = TRUE)
))
predictions_tag_type_pair_ <-
predictions_tag_type_pair %>%
as.data.frame() %>%
mutate(response.level = factor(response.level, levels = c("brood", "fledge", "fail", "aband")))
# Calculate the sample sizes
sample_sizes_ <- stats_breed_pair_brm$data %>%
group_by(tag_type_combo, state) %>%
summarise(n = n()) %>%
ungroup() %>%
rename(x = tag_type_combo,
response.level = state) %>%
bind_rows(., data.frame(x = c("GPS-GPS", "GPS-GPS",
"GPS-control", "GPS-control", "GPS-control",
"PTT-PTT"),
response.level = c("fledge", "aband",
"fledge", "aband", "brood",
"aband"),
n = rep(0, 6)))
# Merge the sample sizes with your predictions
predictions_tag_type_pair_ <- predictions_tag_type_pair_ %>%
left_join(sample_sizes_, by = c("x", "response.level"))
custom_colors_ <- brewer.pal(9, "Set1")[c(9, 3, 2)]
tag_type_breeding_fate_plot_pair <-
ggplot(data = predictions_tag_type_pair_) +
geom_errorbar(aes(x = response.level, y = predicted, ymin = conf.low, ymax = conf.high, color = x),
size = 0.3, linetype = "solid", width = 0.1) +
geom_point(aes(x = response.level, y = predicted, fill = x, shape = response.level), size = 4, color = "grey40", alpha = 1) +
geom_text(aes(label = n, y = -0.2, x = response.level),
size = 3, color = "black", vjust = 0) +
theme(legend.position = "none",
text = element_text(family = "Verdana"),
axis.ticks = element_blank(),
axis.title.y = element_text(size = 11, colour = "grey40"),
axis.title.x = element_blank(),
axis.text.y = element_text(size = 10),
axis.text.x = element_text(size = 10, colour = "grey40", angle = 45, hjust = 1),
plot.margin = unit(c(0.2,0.2,0.2,0.2), "cm"),
panel.background = element_rect(fill = 'transparent', color = NA),
plot.background = element_rect(fill = 'transparent', color = NA),
panel.spacing = unit(1, "cm"),
panel.grid.major.y = element_line(size = 0.25, lineend = "round", colour = "grey90"),
panel.grid.minor = element_blank(),
panel.grid.major.x = element_blank()) +
ylab("fate probability\n(± 95% CrI)") +
xlab("breeding attempt fate") +
scale_y_continuous(limits = c(-0.2, 1), breaks = c(0, 0.25, 0.5, 0.75, 1)) +
scale_x_discrete(labels = c("abandoned", "hatched", "failed", "fledged")) +
facet_grid(~ x) +
scale_shape_manual(values = c(21, 22, 23, 24, 25))
tag_type_breeding_fate_plot_pair
```