Badger data for KK - LS - TP area

Author

Virginia Morera-Pujol

Clean data

Clean the badger data and select only the badgers captured in the Kk-Ls-Tp area

Code

# subset quartiles in kk area
kk_quartiles = st_intersection(quart, counties)

# subset badgers in kk area
kk_badgers <- badgers_vacc %>% 
  filter(QUARTILE %in% c(kk_quartiles$QUARTILE))

# clean badger data
badgers_clean <- kk_badgers %>% 
  mutate(CAPTURE_BLOCK_EVENT = paste(QUARTILE, EVENT_NO, sep = "_")) %>% 
  # select necessary variables
  dplyr::select(SETT_ID, BADGER_ID, DATE_CAUGHT = DATE_CAPTURE, SEX, AGE, WEIGHT,
                VACCINATION, CAPTURE_BLOCK_EVENT,
                ECTOPARASITES, BADGER_STATUS,
                ECTOPARASITE_TYPE, BADGER_ACTION, DATE_ENTERED, 
  ) %>% 
  # modify date caught 
  mutate(SEX = recode(SEX, 'F' = "Female", 'FEMALE' = "Female", 
                      'M' = "Male", 'MALE' = "Male", 
                      'X' = NA_character_), 
         # clean and recode Age
         AGE = recode(AGE, 'adult' = "Adult", 'ADULT' = "Adult",'CUB' = "Cub", 
                      'JUVENILE' = "Juvenile", 'OLD' = "Old", 
                      '9' = NA_character_), 
         # clean and recode badger status
         BADGER_STATUS = recode(BADGER_STATUS, 
                                '2014 badger' = NA_character_, 
                                '2nd Chip' = NA_character_, 
                                'clinical suspec' = NA_character_,
                                'Ann Barber' = NA_character_, 
                                'AnnBarber' = NA_character_,
                                'Deceased' = "Removed",
                                'BADGER' = "Badger",
                                'EPI Removed Under Permission' = "Removed", 
                                'Euthanised' = "Removed", 
                                'Euthanised by V' = "Removed", 
                                'Goodger not on' = "Goodger", 
                                'new' = "New",
                                'NEW' = "New",
                                'Not on system' = NA_character_,
                                'Not on System' = NA_character_,
                                'Old thin' = NA_character_,
                                'poor condition' = NA_character_,
                                'other' = NA_character_,
                                'RECAPTURE' = "Recapture", 
                                'Recapture (14 days)' = "Recapture",
                                'Recapture (30 days)' = "Recapture",
                                'Recaptured' = 'Recapture',
                                'removed' = "Removed",
                                'roadside injury' = NA_character_,
                                'Positive test' = "Removed", 
                                'Positive Test' = "Removed", 
                                'Positve Test' = "Removed",
                                'road kill' = "Removed"), 
         # id programme
         # PROGRAMME = "Vaccination", 
         # clean weight variable
         WEIGHT = na_if(WEIGHT, 0), 
         YEAR = year(DATE_CAUGHT), 
         MONTH = month(DATE_CAUGHT, label = T), 
         BADGER_STATUS = recode(BADGER_STATUS, 
                                'Badger' = "New", 
                                'Goodger' = "Recapture")) %>% 
  filter(BADGER_ACTION == "RELEASED", 
         BADGER_STATUS %in% c("New", "Recapture")) %>%
  droplevels()

Display summaries of the data

Number of new captures and recaptures per year, and recapture percentage

Code

badgers_clean %>% 
  group_by(YEAR, BADGER_STATUS) %>% 
  summarise(N_badgers = n()) %>% 
  ungroup() %>% 
  pivot_wider(names_from = "BADGER_STATUS", values_from = N_badgers) %>% 
  rowwise() %>% 
  mutate(Total = sum(New, Recapture, na.rm = T), 
         Percent_recapture = round(100*(Recapture/Total), 1)) %>% 
  kbl(format = "markdown")

YEAR	New	Recapture	Total	Percent_recapture
2018	78	NA	78	NA
2019	517	248	765	32.4
2020	842	578	1420	40.7
2021	1002	842	1844	45.7
2022	767	783	1550	50.5
2023	773	610	1383	44.1
NA	30	1	31	3.2

Total number of captures and recaptures

Code

ggplot(badgers_clean) + 
  geom_bar(aes(x = BADGER_STATUS, fill = BADGER_ACTION)) + 
  theme_bw() +
  labs(x = "Badger status", fill = "Action", title = "Total badgers captured")

Number of recaptures per year

Code

recap_count <- badgers_clean %>% 
  group_by(BADGER_ID) %>% 
  summarise(n_captures = n(), 
            SETT_ID = head(SETT_ID, 1)) %>% 
  ungroup()

ggplot(recap_count) + 
  geom_bar(aes(x = n_captures)) +
  scale_x_continuous(breaks = 1:10, labels = 1:10) +
  theme_bw() + 
  labs(x = "Number of recaptures", title = "How many times have animals been recaptured?")

Yearly captures and recaptures

Code

ggplot(badgers_clean) + 
  geom_bar(aes(x = YEAR, fill = BADGER_STATUS)) + 
  scale_x_continuous(breaks = 2018:2023, labels = 2018:2023) +
  theme_bw() +
  labs(fill = "Badger status", title = "Badgers captured yearly")

Monthly captures and recaptures

Code

ggplot(badgers_clean) + 
  geom_bar(aes(x = MONTH, fill = BADGER_STATUS)) + 
  # scale_x_continuous(breaks = 2018:2023, labels = 2018:2023) +
  theme_bw() +  
  labs(fill = "Badger status", title = "Badgers captured monthly")

Merge with sett data so we can have spatial information and plot

The badger dataset on its own does not have information on where the badger was captured, but it does have information on what sett they were captured in (or nearby). Therefore, we can combine the information on the sett dataset and obtain coordinates for each badger (that will be sett coordinates)

Overall

Code

badgers_clean_sf <- badgers_clean %>% 
  left_join(sett_all) %>% 
  st_set_geometry("geometry")

ggplot(badgers_clean_sf) + 
  geom_sf(data = counties) + 
  geom_sf(aes(col = factor(YEAR))) + 
  geom_sf(data = ded2, col = "red", fill = NA) + 
  scale_colour_viridis_d() + 
  labs(col = "Year") + 
  theme_bw() +
  labs(col = "Capture year", title = "Badgers captured overall")

Yearly

Code

ggplot(badgers_clean_sf %>% filter(!is.na(YEAR))) + 
  geom_sf(data = counties) + 
  geom_sf(aes(col = BADGER_STATUS)) + 
  geom_sf(data = ded2, fill = NA, col = "red") + 
  scale_colour_viridis_d() + 
  facet_wrap(~YEAR) + 
  theme_bw()  + 
  theme(legend.position="top") + 
  labs(col = "Badger status", title = "Badgers captured yearly")

Effort

The effort data comes from a dataset obtained from DAFM detailing the capture events conducted in each quartile, and the date each commenced and finished, which allows us to know how many days each capture event lasted.

With this dataset we can know how many days of effort were conducted in each quartile. Therefore, the smallest unit of effort (spatially) that we can obtain from this dataset is the quartile.

Effort overall

Code

vac_cap_clean <- vac_cap %>% 
  filter(VACCINE_STATUS == "APPROVED") %>% 
  mutate(DATE_COMMENCED = ymd_hm(DATE_COMMENCED), 
         DATE_COMPLETED = ymd_hm(DATE_COMPLETED),
         DURATION = DATE_COMPLETED- DATE_COMMENCED, 
         YEAR = year(DATE_COMMENCED)) %>% 
  select(QUARTILE, EVENT, CLUSTER_ID, DATE_COMMENCED, DATE_COMPLETED, CHNG_NO, 
         DURATION, YEAR) %>% 
  drop_na(DATE_COMMENCED) %>% 
  distinct() 

vac_cap_sum <- vac_cap_clean %>% 
  mutate(CAPTURE_BLOCK_EVENT_vacc = paste(EVENT, DATE_COMMENCED)) %>% 
  group_by(QUARTILE) %>% 
  summarise(NDAYS = as.numeric(sum(DURATION))) %>% 
  ungroup()

quart_effort <- left_join(kk_quartiles, vac_cap_sum)

ggplot() + 
  geom_sf(data = counties, fill = "lightgray", col = "black") + 
  geom_sf(data = quart_effort, aes(fill = NDAYS)) + 
  geom_sf(data = ded2, fill = NA, col = "red") + 
  # geom_sf(data = sett_all %>% filter(YEAR > 2018), size = 0.5, col = "red") +
  scale_fill_viridis_c(na.value = NA) + 
  theme_bw() + 
  labs(title = "Sampling effort", 
       fill = "N days")

Effort yearly

Code

vac_cap_sum <- vac_cap_clean %>% 
  mutate(CAPTURE_BLOCK_EVENT_vacc = paste(EVENT, DATE_COMMENCED)) %>% 
  group_by(YEAR, QUARTILE) %>% 
  summarise(NDAYS = as.numeric(sum(DURATION))) %>% 
  ungroup() %>% 
  pivot_wider(names_from = "YEAR", values_from = "NDAYS")

quart_effort <- left_join(kk_quartiles, vac_cap_sum) %>% 
  pivot_longer(c("2018", "2019", "2020", "2021", "2022", "2023"), 
               names_to = "YEAR", values_to = "NDAYS")

ggplot() + 
  geom_sf(data = counties, fill = "lightgray", col = "black") + 
  geom_sf(data = quart_effort, aes(fill = NDAYS), col = NA) + 
  # geom_sf(data = sett_all %>% filter(YEAR > 2018), size = 0.5, col = "red") +
  geom_sf(data = ded2, fill = NA, col = "red") + 
  scale_fill_viridis_c(na.value = NA) + 
  theme_bw() + 
  facet_wrap(~YEAR) + 
  labs(title = "Yearly sampling effort", 
       fill = "N days")