In pratikunterwegs/holeybirds: Link bird moult with movement.

Environmental covariates at residence patches

library(data.table)

library(sf)
library(terra)

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

# for gams
library(mgcv)

Extract environmental covariates at patches

data <- fread("data/results/data_patch_summary_ppa.csv")
points <- fread("data/results/data_patch_points_ppa.csv")

points <- points[, c("id", "patch", "x", "y", "date")]
points <- st_as_sf(
  points,
  coords = c("x", "y"),
  crs = 2039
)

# load rasters
ndvi <- rast("data/rasters/raster_hula_ndvi_2039.tif")
vis <- rast("data/rasters/raster_hula_visibility.tif")

ndvi_samp <- terra::extract(
  ndvi,
  vect(points)
)

vis_samp <- terra::extract(
  vis,
  vect(points)
)

points <- st_drop_geometry(points)
setDT(points)

points[, c("ndvi", "vis") := list(
  ndvi_samp$raster_hula_ndvi_2039,
  vis_samp$raster_hula_visibility
)]

rm(vis, ndvi)
gc()

# mean ndvi and vis per patch and date
points <- points[, lapply(.SD, mean), by = c("id", "patch", "date")]

# save
fwrite(
  points,
  "data/results/data_patch_env.csv"
)

Link patch covariates to RRV

# load data patch summary
data <- fread("data/results/data_patch_summary_ppa.csv")
patch_env <- fread("data/results/data_patch_env.csv")

data <- merge(
  data,
  patch_env
)

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

# change name of TAG to id in RRV
setnames(rrv, "TAG", "id")
data <-
  merge(
    data,
    rrv,
    by = c("id", "date"),
    all.x = T
  )

# filter
data <- data[!is.na(sp) & !is.na(treat)]

Duration between patches

data[, duration_bw_patches := (
  (time_start - data.table::shift(time_end)) / 60 # in minutes
), by = c("id", "date")]

# plot histogram of duration between patches
p <- ggplot(data) +
  geom_histogram(
    aes(duration_bw_patches, y = ..density..)
  ) +
  facet_grid(sp ~ treat) +
  labs(
    x = "duration between patches (mins)"
  )
p +
  scale_x_log10()

ggsave(
  p,
  filename = "figures/fig_duration_bw_patches.png"
)

Checking for tiSSA requirement

p <- ggplot(data) +
  geom_jitter(
    aes(
      duration_bw_patches, dist_bw_patch,
      col = treat
    )
  ) +
  geom_hline(
    yintercept = 100,
    lty = 2
  ) +
  geom_vline(
    xintercept = 30,
    lty = 2
  ) +
  scale_x_sqrt() +
  scale_y_sqrt() +
  facet_grid(~sp)

ggsave(
  p,
  filename = "figures/fig_distance_duration_bw_patches.png"
)

Movements between patches

# patch switches
patch_switches <- data[, list(
  N = .N,
  total_dist = sum(dist_bw_patch, na.rm = T),
  mean_duration = mean(duration) / 3600, # in hours
  tracking_duration = (max(time_end) - min(time_start)) / 3600 # in hours
), by = c("id", "date", "sp", "rrv_calc", "treat")]

patch_switches <- patch_switches[!is.na(sp)]

# convert to patch switches per hour
patch_switches[, dist_ph := total_dist / tracking_duration]

# remove NA values
patch_switches <- patch_switches[!is.na(rrv_calc)]

# asign sp as factor
patch_switches[, sp := factor(sp, levels = c(
  "Pycnonotus", "Passer", "Acrocephalus"
))]

Examine the number of patch switches in relation to RRV.

ggplot(patch_switches) +
  geom_jitter(
    aes(
      rrv_calc, N / tracking_duration,
      col = treat
    )
  ) +
  geom_smooth(
    aes(
      rrv_calc, N / tracking_duration
    ),
    method = "gam",
    formula = as.formula(
      "y ~ s(x, k = 3)"
    )
  ) +
  facet_grid(~sp)

Model patch switches by RRV

# scale patch switches by tracking duration
patch_switches[, n_per_hr := N / tracking_duration]

# fit gam
mod_ps <- gam(
  n_per_hr ~ s(rrv_calc, k = 3, by = sp) +
    s(sp, bs = "re"),
  data = patch_switches
)

# visualise gam
gratia::draw(mod_ps)

summary(mod_ps)

# write summary to file
writeLines(
  capture.output(
    summary(mod_ps)
  ),
  con = "data/results/mod_summary_rrv_n_patches.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(mod_ps),
  custom.model.names = c("Coefficients"),
  single.row = T,
  dcolumn = T,
  custom.coef.names = c(
    "Intercept",
    "WGI - bulbul", "WGI - sparrow", "WGI - warbler",
    "Species"
  ),
  caption = "Generalised additive model coefficients for residence patch switches."
)

# for word file
texreg::wordreg(
  l = list(mod_ps),
  file = "docs/model_patch_switches.docx",
  custom.model.names = c("Coefficients"),
  single.row = T,
  dcolumn = T,
  custom.coef.names = c(
    "Intercept",
    "WGI - bulbul", "WGI - sparrow", "WGI - warbler",
    "Species"
  ),
  caption = "Generalised additive model coefficients for residence patch switches."
)

# make prediction table
pred_data_patches_rrv <- CJ(
  sp = as.factor(unique(data$sp)),
  id = "new",
  treat = as.factor(unique(data$treat)),
  rrv_calc = seq(0, 20, 0.5)
)

# filter for unrealistic
pred_data_patches_rrv <- pred_data_patches_rrv[(
  (treat == "NonMoulting" & rrv_calc == 0) |
    (treat == "Moulting" & between(rrv_calc, 2, 12)) |
    (treat == "Manipulated" & between(rrv_calc, 12, 20))
)]

# get prediction
pred_ps <- predict(
  mod_ps,
  newdata = pred_data_patches_rrv,
  allow.new.levels = T, se.fit = T
)

pred_data_patches_rrv$pred <- pred_ps$fit
pred_data_patches_rrv$se <- pred_ps$se.fit

Plot patch switches per hour

# prepare summary statistics for switches per hour and wing gap index
patch_switches_summary <- patch_switches[, unlist(
  lapply(.SD, function(x) {
    list(
      mean = mean(x),
      sd = sd(x)
    )
  }),
  recursive = FALSE
), .SDcols = c("rrv_calc", "n_per_hr"), by = c("sp", "treat")]

Plotting code for patch switches per hour of tracking, in relation to wing gap index.

# save as object and write intermediate to file
fig_patch_switches <-
  ggplot(patch_switches) +
  geom_jitter(
    aes(
      rrv_calc, n_per_hr,
      col = treat
    ),
    shape = 1
  ) +
  geom_ribbon(
    data = pred_data_patches_rrv[!(sp == "Pycnonotus" & rrv_calc > 17) &
      !sp == "Acrocephalus"],
    aes(
      rrv_calc,
      ymin = pred - se,
      ymax = pred + se
    ),
    fill = "transparent",
    col = "grey",
    lty = 1
  ) +
  geom_line(
    data = pred_data_patches_rrv[!(sp == "Pycnonotus" & rrv_calc > 17) &
      !sp == "Acrocephalus"],
    aes(
      rrv_calc, pred
    ),
    col = "indianred"
  ) +
  geom_linerange(
    data = patch_switches_summary,
    aes(
      rrv_calc.mean,
      ymin = n_per_hr.mean - n_per_hr.sd,
      ymax = n_per_hr.mean + n_per_hr.sd
    )
  ) +
  geom_linerange(
    data = patch_switches_summary,
    aes(
      rrv_calc.mean, n_per_hr.mean,
      xmin = rrv_calc.mean - rrv_calc.sd,
      xmax = rrv_calc.mean + rrv_calc.sd
    )
  ) +
  geom_point(
    data = patch_switches_summary,
    aes(
      rrv_calc.mean, n_per_hr.mean,
      fill = treat
    ),
    shape = 21,
    size = 2
  ) +
  facet_grid(
    ~sp
  ) +
  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"
  ) +
  facet_grid(
    cols = vars(sp),
    scales = "free",
    labeller = labeller(
      sp = c(
        "Acrocephalus" = "Clamorous\nreed warbler",
        "Passer" = "House sparrow",
        "Pycnonotus" = "White-spectacled\nbulbul"
      )
    )
  ) +
  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 →)",
    # x = NULL,
    y = "Patch switches / hr",
    fill = NULL
  )

# save to file
ggsave(
  fig_patch_switches,
  filename = "figures/fig_patch_switches.png",
  width = 120,
  height = 70,
  units = "mm"
)

Model inter-patch distance per hour

# plot rrv with errors per treatment, and patch distance w/ errors per trt
psdf <- patch_switches[, unlist(lapply(.SD, function(x) {
  list(
    mean = as.numeric(mean(x)),
    sd = sd(x, na.rm = T)
  )
}), recursive = F),
.SDcols = c("dist_ph", "rrv_calc"), by = c("sp", "treat")
]

patch_switches[, c("sp", "treat") := lapply(
  .SD, as.factor
), .SDcols = c("sp", "treat")]

mod1 <- gam(
  dist_ph ~ s(rrv_calc, by = sp, k = 3) +
    s(sp, bs = "re"),
  # s(sp, rrv_calc, bs = "re"),
  data = patch_switches
)

summary(mod1)
gratia::draw(mod1)

# save model summary
mod_summary <- summary(mod1)
writeLines(
  capture.output(
    mod_summary
  ),
  con = "data/results/mod_summary_rrv_movement.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",
    "WGI - bulbul", "WGI - sparrow", "WGI - warbler",
    "Species"
  ),
  caption = "Generalised additive model coefficients for distance between residence patches."
)

# for word file
texreg::wordreg(
  list(mod1),
  file = "docs/model_patch_distances.docx",
  custom.model.names = c("Coefficients"),
  single.row = T,
  dcolumn = T,
  custom.coef.names = c(
    "Intercept",
    "WGI - bulbul", "WGI - sparrow", "WGI - warbler",
    "Species"
  ),
  caption = "Generalised additive model coefficients for distance between residence patches."
)

Summarise between-patch movements

# summarise patch switches by species and treatments
patch_distance_summary <- patch_switches[treat != "Manipulated" |
  (treat == "Manipulated" & rrv_calc > 12), list(
  mean_dist_ph = mean(dist_ph, na.rm = TRUE),
  sd_dist_ph = sd(dist_ph, na.rm = TRUE)
), by = c("sp", "treat")]

# save small data
fwrite(
  patch_distance_summary,
  file = "data/results/data_patch_distance_summary.csv"
)

Plot movement per hour over wing gap index

# make prediction table
pred_data <- CJ(
  sp = as.factor(unique(data$sp)),
  id = "new",
  treat = as.factor(unique(data$treat)),
  rrv_calc = seq(0, 20, 0.5)
)

# filter for unrealistic
pred_data <- pred_data[(
  (treat == "NonMoulting" & rrv_calc == 0) |
    (treat == "Moulting" & between(rrv_calc, 2, 12)) |
    (treat == "Manipulated" & between(rrv_calc, 12, 20))
)]

# 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(patch_switches[dist_ph < 900]) +
  # geom_jitter(
  #   aes(
  #     rrv_calc, dist_ph,
  #     col = treat
  #   ),
  #   shape = 1
  # )+
  geom_point(
    data = pred_data,
    aes(
      rrv_calc, pred
    )
  ) +
  facet_grid(
    ~sp
  )

# set factor order of species names
psdf[, sp := factor(sp, levels = c(
  "Pycnonotus", "Passer", "Acrocephalus"
))]

# plot figure
fig_moult_move <-
  ggplot(psdf) +
  geom_hline(
    yintercept = 100,
    col = "grey",
    lty = 2
  ) +
  geom_jitter(
    data = patch_switches[dist_ph < 900],
    aes(
      rrv_calc, dist_ph,
      col = treat
    ),
    shape = 1
  ) +
  geom_ribbon(
    data = pred_data[!(sp == "Pycnonotus" & rrv_calc > 17) &
      !sp == "Acrocephalus"],
    aes(
      rrv_calc,
      ymin = pred - se,
      ymax = pred + se
    ),
    fill = "transparent",
    col = "grey",
    lty = 1
  ) +
  geom_line(
    data = pred_data[!(sp == "Pycnonotus" & rrv_calc > 17) &
      !sp == "Acrocephalus"],
    aes(
      rrv_calc, pred
    ),
    col = "indianred"
  ) +
  geom_linerange(
    aes(
      rrv_calc.mean,
      ymin = dist_ph.mean - dist_ph.sd,
      ymax = dist_ph.mean + dist_ph.sd
    )
  ) +
  geom_linerange(
    aes(
      rrv_calc.mean, dist_ph.mean,
      xmin = rrv_calc.mean - rrv_calc.sd,
      xmax = rrv_calc.mean + rrv_calc.sd
    )
  ) +
  geom_point(
    aes(
      rrv_calc.mean, dist_ph.mean,
      fill = treat
    ),
    shape = 21,
    size = 2
  ) +
  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_y_sqrt()+
  facet_grid(
    cols = vars(sp),
    scales = "free",
    labeller = labeller(
      sp = c(
        "Acrocephalus" = "Clamorous\nreed warbler",
        "Passer" = "House sparrow",
        "Pycnonotus" = "White-spectacled\nbulbul"
      )
    )
  ) +
  coord_cartesian(
    ylim = c(0, NA)
  ) +
  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 = "Between-patch distance (m)",
    fill = NULL
  )

ggsave(
  fig_moult_move,
  filename = "figures/fig_moult_move.png",
  width = 120,
  height = 70,
  units = "mm"
)

Patch durations in relation to wing gap index

# explore the data
ggplot(data) +
  geom_jitter(
    aes(
      vis, duration / (3600),
      col = treat
    )
  ) +
  facet_grid(treat ~ sp) +
  coord_cartesian(
    ylim = c(0, 2)
  )

# handle variable as factor
data[, sp := factor(sp, levels = c(
  "Pycnonotus", "Passer", "Acrocephalus"
))]

# convert to hours
data[, duration_hr := duration / (3600)]

# fit a gam with few knots
mod_duration_rrv <- gam(
  duration_hr ~ s(rrv_calc, by = sp, k = 3) +
    vis + ndvi +
    s(sp, bs = "re"),
  data = data
)

# plot model
gratia::draw(mod_duration_rrv)

# examine model summary
summary(mod_duration_rrv)

# save model summary
mod_summary_duration <- summary(mod_duration_rrv)
writeLines(
  capture.output(
    mod_summary_duration
  ),
  con = "data/results/mod_summary_duration_rrv.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(mod_duration_rrv),
  custom.model.names = c("Coefficients"),
  single.row = T,
  dcolumn = T,
  custom.coef.names = c(
    "Intercept", "Visibility index", "NDVI",
    "WGI - bulbul", "WGI - sparrow", "WGI - warbler",
    "Species"
  ),
  caption = "Generalised additive model coefficients for residence patch duration."
)

# for word file
texreg::wordreg(
  list(mod_duration_rrv),
  file = "docs/model_patch_duration.docx",
  custom.model.names = c("Coefficients"),
  single.row = T,
  dcolumn = T,
  custom.coef.names = c(
    "Intercept", "Visibility index", "NDVI",
    "WGI - bulbul", "WGI - sparrow", "WGI - warbler",
    "Species"
  ),
  caption = "Generalised additive model coefficients for residence patch duration."
)

Plot durations in relation to wing gap index

# prepare predictor data
# make prediction table
pred_data_duration_rrv <- CJ(
  sp = as.factor(unique(data$sp)),
  id = "new",
  treat = as.factor(unique(data$treat)),
  rrv_calc = seq(0, 20, 0.5),
  ndvi = 0.4, # taken from the data mean
  vis = 0.33 # from the data mean
)

# filter for unrealistic
pred_data_duration_rrv <- pred_data_duration_rrv[(
  (treat == "NonMoulting" & rrv_calc == 0) |
    (treat == "Moulting" & between(rrv_calc, 2, 12)) |
    (treat == "Manipulated" & between(rrv_calc, 12, 20))
)]

# get prediction
pred_duration <- predict(
  mod_duration_rrv,
  newdata = pred_data_duration_rrv,
  allow.new.levels = T, se.fit = T
)

pred_data_duration_rrv$pred <- pred_duration$fit
pred_data_duration_rrv$se <- pred_duration$se.fit

# make summary table
patch_duration_summary <- data[, unlist(
  lapply(.SD, function(x) {
    list(
      mean = mean(x, na.rm = T),
      sd = sd(x, na.rm = T)
    )
  }),
  recursive = F
), .SDcols = c("duration_hr", "rrv_calc"), by = c("sp", "treat"), ]

# prepare figure
fig_patch_duration_rrv <-
  ggplot(data) +
  geom_jitter(
    aes(
      rrv_calc,
      duration_hr,
      col = treat
    ),
    shape = 1,
    size = 0.5
  ) +
  geom_ribbon(
    data = pred_data_duration_rrv[!(sp == "Pycnonotus" & rrv_calc > 17) &
      !sp == "Passer"],
    aes(
      rrv_calc,
      ymin = pred - se,
      ymax = pred + se
    ),
    fill = "transparent",
    col = "grey",
    lty = 1
  ) +
  geom_line(
    data = pred_data_duration_rrv[!(sp == "Pycnonotus" & rrv_calc > 17) &
      !sp == "Passer"],
    aes(
      rrv_calc, pred
    ),
    col = "indianred"
  ) +
  geom_linerange(
    data = patch_duration_summary,
    aes(
      rrv_calc.mean,
      ymin = duration_hr.mean - duration_hr.sd,
      ymax = duration_hr.mean + duration_hr.sd
    )
  ) +
  geom_linerange(
    data = patch_duration_summary,
    aes(
      rrv_calc.mean, duration_hr.mean,
      xmin = rrv_calc.mean - rrv_calc.sd,
      xmax = rrv_calc.mean + rrv_calc.sd
    )
  ) +
  geom_point(
    data = patch_duration_summary,
    aes(
      rrv_calc.mean, duration_hr.mean,
      fill = treat
    ),
    shape = 21,
    size = 2
  ) +
  coord_cartesian(
    ylim = c(0, 5)
  ) +
  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_y_sqrt()+
  facet_grid(
    cols = vars(sp),
    scales = "free",
    labeller = labeller(
      sp = c(
        "Acrocephalus" = "Clamorous\nreed warbler",
        "Passer" = "House sparrow",
        "Pycnonotus" = "White-spectacled\nbulbul"
      )
    )
  ) +
  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 = "Duration in patch (hrs)",
    fill = NULL
  )

ggsave(
  fig_patch_duration_rrv,
  filename = "figures/fig_patch_duration_rrv.png",
  width = 120, height = 70, units = "mm"
)

pratikunterwegs/holeybirds documentation built on Aug. 6, 2023, 8:53 a.m.