In pratikunterwegs/holeybirds: Link bird moult with movement.

Selection for sheltered areas by bulbuls, sparrows, and reed warblers

library(readr)
library(dplyr)
library(tidyr)
library(glue)

library(survival)

library(ggplot2)
library(ggtext)
library(colorspace)
library(patchwork)

library(mgcv)

# read in saved data
data_step_covs <- read_csv("data/results/data_ssf_step_covs.csv")

# remove cols
data_step_covs <- select(
  data_step_covs,
  id,
  date,
  matches(c("x", "y")),
  matches(c("ndvi", "vis")),
  case_, step_id_,
  B, O, R, `T`, W,
  -`...1`
)

# read rrv
rrv <- read_csv("data/results/data_daily_rrv.csv")
rrv$date <- as.character(rrv$date)

data_step_covs <- mutate(
  data_step_covs,
  date = as.character(date)
)

data_step_covs <- left_join(
  data_step_covs,
  rrv,
  by = c("id" = "TAG", "date")
)

# long term tracking
df_cov_lt <- filter(
  data_step_covs,
  is.na(rrv_calc),
  sp == "Passer"
)

Summarise by case and treatment

# filter for high ndvi values
data_sel <- filter(
  data_step_covs,
  between(ndvi_end, 0, 1)
)

# filter out data without species or rrv
data_sel <- filter(
  data_sel,
  !is.na(sp),
  !is.na(treat),
  !is.na(rrv_calc)
)

Summarise number of fixes in real and alternate patches.

data_sel %>%
  group_by(
    id, treat,
    case_
  ) %>%
  count() %>%
  group_by(
    case_
  ) %>%
  summarise(
    across(n, .fns = list(mean = mean, sd = sd))
  )

Prepare data for conditional logistic regression.

# select useful columns
sdf <- filter(
  data_sel
) |>
  select(
    sp, treat, case_,
    rrv_calc,
    matches(c("vis", "ndvi"))
  )

sdf <- sdf |>
  group_by(
    sp, treat, case_
  ) |>
  summarise(
    across(
      matches(c("vis", "ndvi", "rrv")),
      .fns = list(
        mean = function(x) mean(x, na.rm = T),
        sd = function(x) sd(x, na.rm = T)
      )
    ),
    .groups = "keep"
  ) |>
  rename(
    vis_mean = vis_end_mean,
    vis_sd = vis_end_sd,
    ndvi_mean = ndvi_end_mean,
    ndvi_sd = ndvi_end_sd
  ) |>
  mutate(
    case_ = if_else(case_, "real", "alternate")
  )

# save
write_csv(
  sdf,
  file = "data/results/data_visibility_steps_summary.csv"
)

Model visibility with RRV

# convert sp and treat to factor
data_sel <- mutate(
  data_sel,
  across(
    c(sp, treat, id),
    .fns = as.factor
  )
)

# fit gam
mod1 <- gam(
  vis_end ~ s(rrv_calc, by = sp, k = 3) +
    s(ndvi_end, k = 5) +
    s(sp, bs = "re"),
  data = filter(
    data_sel, case_
  )
)

summary(mod1)
gratia::draw(mod1)

# save model summary
mod_summary <- summary(mod1)
writeLines(
  capture.output(
    mod_summary
  ),
  con = "data/results/mod_summary_rrv_visibility.txt"
)

This code below is only to create clean outputs for the LaTeX supplementary material.

# get model as tex code for supplement
texreg::texreg(
  list(mod1),
  custom.model.names = c("Coefficients"),
  single.row = T,
  dcolumn = T,
  custom.coef.names = c(
    "Intercept", "NDVI",
    "WGI - bulbul", "WGI - sparrow", "WGI - warbler",
    "Species"
  ),
  caption = "Generalised additive model coefficients for residence patch visibility."
)

texreg::wordreg(
  list(mod1),
  file = "docs/model_patch_vis.docx",
  custom.model.names = c("Coefficients"),
  single.row = T,
  dcolumn = T,
  custom.coef.names = c(
    "Intercept", "NDVI",
    "WGI - bulbul", "WGI - sparrow", "WGI - warbler",
    "Species"
  ),
  caption = "Generalised additive model coefficients for residence patch visibility."
)

Model visibility across RRV

# make prediction table
pred_data <- crossing(
  sp = as.factor(unique(data_sel$sp)),
  id = "new",
  ndvi_end = 0.3,
  rrv_calc = seq(0, 20, 0.5)
)

# get prediction
pred <- predict(mod1, newdata = pred_data, allow.new.levels = T, se.fit = T)

pred_data$pred <- pred$fit
pred_data$se <- pred$se.fit

# explore
ggplot(data_sel) +
  # geom_jitter(
  #   aes(
  #     rrv_calc, vis_mean_end,
  #     col = treat
  #   ),
  #   shape = 1
  # )+
  geom_point(
    data = pred_data,
    aes(
      rrv_calc, pred
    )
  ) +
  facet_grid(
    ~sp
  )

Step selection based on visibility

# split the data by species, first making a copy
data_clogit <- data_sel

# fix factor levels
data_clogit <- mutate(
  data_clogit,
  treat = factor(
    treat,
    levels = c("NonMoulting", "Moulting", "Manipulated")
  ),
  vis_end_inv = 1 - vis_end
)

# now nest the data
data_clogit <- nest(
  data_clogit,
  data = !c("sp", "treat")
)

# apply a gam over the nested data
data_clogit <- mutate(
  data_clogit,
  mod_fit = lapply(
    data, function(df) {
      survival::clogit(
        formula = case_ ~ (vis_end) + ndvi_end +
          strata(step_id_),
        method = "approximate",
        data = df
      )
    }
  )
)

Extract model coefficients for supplementary material.

texreg::texreg(
  l = data_clogit$mod_fit,
  dcolumn = T,
  booktabs = T,
  custom.header = list(
    "White-spectacled bulbul" = 1:3,
    "House sparrow" = 4:6,
    "Clamorous reed warbler" = 7:9
  ),
  custom.model.names = rep(
    c("Non-moulting", "Molting", "Manipulated"),
    times = 3
  ),
  custom.coef.names = c("Visibility index", "NDVI"),
  caption = "Log relative selection strengths (RSS) for patch visibility and NDVI."
)

texreg::wordreg(
  data_clogit$mod_fit,
  file = "docs/model_ssf.docx",
  dcolumn = T,
  booktabs = T,
  custom.header = list(
    "White-spectacled bulbul" = 1:3,
    "House sparrow" = 4:6,
    "Clamorous reed warbler" = 7:9
  ),
  custom.model.names = rep(
    c("Non-moulting", "Molting", "Manipulated"),
    times = 3
  ),
  custom.coef.names = c("Visibility index", "NDVI"),
  caption = "Log relative selection strengths (RSS) for patch visibility and NDVI."
)

# get model outputs
data_clogit <- mutate(
  data_clogit,
  mod_fit = lapply(mod_fit, broom::tidy)
)

# unnest model fits
data_clogit <- unnest(
  data_clogit,
  cols = "mod_fit"
)

# mark sig
data_clogit <- mutate(
  data_clogit,
  sig = p.value < 0.05
)

# remove data col and save
data_clogit <- select(data_clogit, -data)
write_csv(
  data_clogit,
  file = "data/results/data_clogit_fit_vis.csv"
)

# merge rrv mean data and clogit fits
data_clogit <- left_join(
  data_clogit,
  filter(sdf, case_ == "real")
)

# add text label
data_clogit <- mutate(
  data_clogit,
  p_val_round = round(p.value, digits = 3),
  label = glue(
    "p = {p_val_round}"
  ),
  # handle if p is very small
  label = if_else(
    p_val_round < 0.0001,
    glue("p < 0.001"),
    label
  ),
  # add the z statistic
  label = glue(
    "ß = {round(estimate, 2)}
    {label}"
  )
)

# arrange data by bird name
data_clogit <- arrange(
  data_clogit,
  desc(sp)
)

Plot comparison of real and potential patches visibility

# set factor levels
sdf <- mutate(
  sdf,
  sp = factor(sp, levels = c(
    "Pycnonotus", "Passer", "Acrocephalus"
  ))
)

# set factors
data_clogit <- mutate(
  data_clogit,
  sp = factor(sp, levels = c(
    "Pycnonotus", "Passer", "Acrocephalus"
  ))
)

# make figure selection coefficients
fig_vis_selection <-
  ggplot(sdf) +
  geom_jitter(
    data = filter(
      data_sel,
      !case_
    ),
    aes(
      rrv_calc, vis_end
    ),
    col = "lightgrey",
    shape = 3,
    size = 0.2
  ) +
  geom_jitter(
    data = filter(
      data_sel,
      case_
    ),
    aes(
      rrv_calc, vis_end,
      col = treat,
    ),
    size = 1,
    shape = 1
  ) +
  geom_linerange(
    aes(
      x = rrv_calc_mean,
      ymin = vis_mean - vis_sd,
      ymax = vis_mean + vis_sd
    ),
    position = position_nudge(
      x = c(-1, 1)
    )
  ) +
  geom_linerange(
    aes(
      x = rrv_calc_mean, y = vis_mean,
      xmin = rrv_calc_mean - rrv_calc_sd,
      xmax = rrv_calc_mean + rrv_calc_sd
    ),
    position = position_nudge(
      x = c(-1, 1)
    )
  ) +
  geom_point(
    aes(
      rrv_calc_mean, vis_mean,
      shape = case_,
      fill = treat
    ),
    size = 3,
    position = position_nudge(
      x = c(-1, 1)
    )
  ) +
  geom_text(
    data = filter(
      data_clogit,
      term == "vis_end"
    ) %>% mutate(
      x_pos = case_when(
        treat == "NonMoulting" ~ 0,
        treat == "Moulting" ~ 6,
        treat == "Manipulated" ~ 12,
        T ~ NA_real_
      )
    ),
    aes(
      x = x_pos, # manual placement
      y = 1.1, # manual placement,
      label = label,
      col = treat
    ),
    hjust = 0,
    fontface = "italic",
    size = 3
  ) +
  facet_grid(
    cols = vars(sp),
    scales = "free_x",
    labeller = labeller(
      sp = c(
        "Acrocephalus" = "Clamorous\nreed warbler",
        "Passer" = "House sparrow",
        "Pycnonotus" = "White-spectacled\nbulbul"
      )
    )
  ) +
  scale_fill_discrete_sequential(
    palette = "Batlow",
    l1 = 30, l2 = 60,
    breaks = c("NonMoulting", "Moulting", "Manipulated"),
    labels = c("Non-molting", "Molting", "Manipulated")
  ) +
  scale_colour_discrete_sequential(
    palette = "Batlow",
    l1 = 50, l2 = 50,
    guide = "none"
  ) +
  scale_shape_manual(
    values = c(
      "real" = 21,
      "alternate" = 25
    ),
    limits = c("alternate", "real"),
    # breaks = c("real", "alternate"),
    labels = c("Potential", "Real")
  ) +
  scale_y_continuous(
    breaks = seq(0, 1, 0.25)
  ) +
  theme_test(
    base_size = 10,
    base_family = "Arial"
  ) +
  theme(
    legend.position = "top",
    strip.background = element_blank(),
    strip.text = element_text(
      face = "italic"
    )
  ) +
  labs(
    x = "Wing gap index (More gappy wing →)",
    y = "Visibility of residence patches",
    fill = NULL,
    shape = NULL
  ) +
  coord_cartesian(
    ylim = c(0, 1.15)
  ) +
  guides(
    fill = guide_legend(
      override.aes = list(
        shape = 22
      )
    ),
    shape = guide_legend(
      override.aes = list(
        size = 2.5
      )
    )
  )

ggsave(
  fig_vis_selection,
  filename = "figures/fig_vis_selection.png",
  height = 90,
  width = 180,
  units = "mm"
)

Compare real and potential patches NDVI

fig_ndvi_selection <- ggplot(sdf) +
  geom_jitter(
    data = filter(
      data_sel,
      !case_
    ),
    aes(
      rrv_calc, ndvi_end
    ),
    col = "lightgrey",
    shape = 4,
    size = 0.2
  ) +
  geom_jitter(
    data = filter(
      data_sel,
      case_
    ),
    aes(
      rrv_calc, ndvi_end,
      col = treat,
    ),
    shape = 1
  ) +
  geom_linerange(
    aes(
      x = rrv_calc_mean,
      ymin = ndvi_mean - ndvi_sd,
      ymax = ndvi_mean + ndvi_sd
    ),
    position = position_nudge(
      x = c(-1, 1)
    )
  ) +
  geom_linerange(
    aes(
      x = rrv_calc_mean, y = ndvi_mean,
      xmin = rrv_calc_mean - rrv_calc_sd,
      xmax = rrv_calc_mean + rrv_calc_sd
    ),
    position = position_nudge(
      x = c(-1, 1)
    )
  ) +
  geom_point(
    aes(
      rrv_calc_mean, ndvi_mean,
      shape = case_,
      fill = treat
    ),
    size = 3,
    position = position_nudge(
      x = c(-1, 1)
    )
  ) +
  scale_fill_discrete_sequential(
    palette = "Batlow",
    l1 = 30, l2 = 60,
    breaks = c("NonMoulting", "Moulting", "Manipulated"),
    labels = c("N. Moulting", "Moulting", "Manipulated")
  ) +
  scale_colour_discrete_sequential(
    palette = "Batlow",
    l1 = 50, l2 = 50,
    guide = "none"
  ) +
  scale_shape_manual(
    values = c(
      "real" = 21,
      "alternate" = 25
    ),
    limits = c("alternate", "real"),
    # breaks = c("real", "alternate"),
    labels = c("Potential", "Real")
  ) +
  facet_grid(
    cols = vars(sp),
    scales = "free_x"
  ) +
  theme_test(
    base_size = 10,
    base_family = "Arial"
  ) +
  theme(
    legend.position = "top",
    strip.background = element_blank(),
    strip.text = element_text(
      face = "italic"
    )
  ) +
  labs(
    x = "RRV (More ragged wing)",
    y = "NDVI at next patch",
    fill = NULL,
    shape = NULL
  ) +
  guides(
    fill = guide_legend(
      override.aes = list(
        shape = 22
      )
    ),
    shape = guide_legend(
      override.aes = list(
        size = 3
      )
    )
  )

ggsave(
  fig_ndvi_selection,
  filename = "figures/fig_ndvi_sel_rrv.png",
  height = 87,
  width = 178,
  units = "mm"
)

Patch selection over the long term

df_cov_lt

sdf <- filter(
  df_cov_lt
) |>
  select(
    sp, treat, case_,
    date, id,
    RRV,
    matches(c("vis", "ndvi"))
  )

sdf <- mutate(
  sdf,
  date = as.Date(date)
) |>
  group_by(
    id
  ) |>
  mutate(
    days = as.numeric(date - min(date))
  )

sdf <- sdf |>
  group_by(
    treat, case_, RRV, days, id
  ) |>
  summarise(
    across(
      matches(c("vis", "ndvi")),
      .fns = list(
        mean = function(x) mean(x, na.rm = T),
        sd = function(x) sd(x, na.rm = T)
      )
    ),
    .groups = "keep"
  ) |>
  rename(
    vis_mean = vis_end_mean,
    vis_sd = vis_end_sd,
    ndvi_mean = ndvi_end_mean,
    ndvi_sd = ndvi_end_sd
  ) |>
  mutate(
    case_ = if_else(case_, "real", "alternate")
  )

ggplot(sdf) +
  stat_summary(
    aes(
      plyr::round_any(days, 3),
      ndvi_mean,
      col = as.factor(RRV),
      shape = case_
    )
  ) +
  scale_color_viridis_d(
    option = "H"
  ) +
  facet_wrap(
    ~id
  )