--- 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 <- function (x, y, x.0 = mean (x)) {<- lm (log (y) ~ log (x))<- rlm (log (y) ~ log (x), method = "M" )<- coef (sma (log (y) ~ log (x)))[2 ]<- coef (sma (log (y) ~ log (x), robust = T))[2 ]<- y * (x.0 / x) ^ b.msa.ols<- y * (x.0 / x) ^ b.msa.rob<- data.frame (SMI.ols, SMI.rob, x, y)<- data.table (x = seq (min (x), max (x), length = 100 )) %>% : = predict (logM.ols, newdata = .) %>% exp] %>% : = 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} = 0.29 = 0.08 = 0.16 = 1.8 = 2.0 = 1.2 ``` ## PCA of structural traits ```{r} <- read.csv (here ("data/body_traits_season.csv" ))``` #### check most relevent structural traits with PCA ```{r} <- %>% group_by (ring) %>% arrange (measure) %>% slice (1 ) %>% ungroup () %>% :: 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 <- %>% 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 ```{r} <- %>% 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 ```{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 ``` [ defined in chunk above ] to derive the SMI for each body mass measure (see Eq 2 of Peig and Green (2009))```{r} <- scaledMassIndex (x = body_traits_avg %>% ungroup () %>% :: select (weight, wing) %>% na.omit () %>% pull (wing), y = body_traits_avg %>% ungroup () %>% :: select (weight, wing) %>% na.omit () %>% pull (weight))<- %>% filter (! is.na (weight)) %>% bind_cols (., smi_) %>% :: select (ring, tag_type, tag_type_, treatment_group, measure, SMI.ols, weight, day_of_year, date_diff) %>% rename (smi_wing = SMI.ols)``` ```{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))%>% filter (! is.na (date_diff)) %>% pull (date_diff) %>% min ()%>% 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} <- lmer (smi_wing ~ tag_type * measure + scale (day_of_year) + (1 | ring), data = body_traits_smi, REML = FALSE )<- 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} <- 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} <- 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} <- 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} <- 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} <- 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 ```{r} %>% 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))``` ```{r} %>% group_by (tag_type_combo) %>% summarise (n_nests = n_distinct (nest_id))``` ```{r} %>% :: 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 ``` 
