Diego J. Lizcano 
Jorge Velásquez-Tibata
Can we use abundance to calculate Bird Friendliness Index (BFI) ?
Load libraries and set some options.
# set up
library(cluster)
library(NbClust)
library(vegan)
library(unmarked)
library(mgcv)
library(mgcViz)
library(tidyverse)
library(gt)
library(lubridate)
library(stringr)
library(readxl)
# setwd("C:/CodigoR/AudubonPanama/R/BFI")
options(scipen=99999)
options(max.print=99999)
options(stringsAsFactors=F)
attributes <- read_excel("C:/CodigoR/AudubonPanama/data/species_data_4_bfi.xlsx", sheet = "attributes")
abundance <- read_excel("C:/CodigoR/AudubonPanama/data/species_data_4_bfi.xlsx", sheet = "abundance_result")We join the different data tables and calculate BFI components per count site.
# diss matirx
print(names(attributes))
#> [1] "sp" "Scientific Name (ACAD taxonomy)"
#> [3] "Scientific Name (Elton traits)" "Group"
#> [5] "CCS_max" "Diet-Inv"
#> [7] "Diet-Vend" "Diet-Vect"
#> [9] "Diet-Vfish" "Diet-Vunk"
#> [11] "Diet-Scav" "Diet-Fruit"
#> [13] "Diet-Nect" "Diet-Seed"
#> [15] "Diet-PlantO" "ForStrat-watbelowsurf"
#> [17] "ForStrat-wataroundsurf" "ForStrat-ground"
#> [19] "ForStrat-understory" "ForStrat-midhigh"
#> [21] "ForStrat-canopy" "ForStrat-aerial"
#> [23] "Nocturnal" "BodyMass-Value"
d_gower <- daisy(attributes[,-c(1,2,3,4,5)], # remove some columns
metric="gower", stand=TRUE,
weights=c(rep(3/9/10, 10),
rep(3/9/8, 8),
rep(3/9/1, 1))) # pesos por atributo son 19
# clusters
clust <- NbClust(diss=d_gower, distance=NULL, method="ward.D",
index="frey") # silhouette before
#>
#> Only frey, mcclain, cindex, sihouette and dunn can be computed. To compute the other indices, data matrix is needed
# best clusters
inddf <- attributes %>% select(sp) %>%
mutate(func_spec=clust$Best.partition) %>%
arrange(func_spec)
# write.csv (inddf %>% left_join(d4 %>% distinct(species, eng_name)),
# "func_species_groupings.csv", row.names=F)
# see table
# gt(inddf) %>%
# tab_header(
# title = md("Funtional group **Parita birds**")#,
# # subtitle = md("`gtcars` is an R dataset")
# ) |> opt_interactive()
#
datatable(
inddf, extensions = 'Buttons', options = list(
dom = 'Bfrtip',
buttons = c('copy', 'csv')
)
)# join functional species and conservation scores to max counts (mean)
score_wt_counts <- abundance %>% left_join(attributes) %>% left_join(inddf) %>%
mutate(count_score=mean*CCS_max) %>%
select(site, count_score) %>%
group_by(site) %>%
summarise(tot_wt_count=sum(count_score)) %>%
ungroup()
# DT::datatable(score_wt_counts)
# # waterbird
# # join functional species and conservation scores to max counts
# score_wt_counts_wb <- abundance %>% left_join(attributes) %>% left_join(inddf) %>%
# filter(Group=="waterbird") %>%
# mutate(count_score=mean*CCS_max) %>%
# # select(farm, pc_station_id, month_id, count_score) %>%
# group_by(site) %>%
# summarise(tot_wt_count_wb=sum(count_score)) %>%
# ungroup()
#
#
# calc shannon, exchange max_count by mean
shan_div <- abundance %>% left_join(attributes) %>% left_join(inddf) %>%
select(site, func_spec, mean) %>%
group_by(site, func_spec) %>%
summarise(tot_count=sum(mean)) %>%
ungroup() %>%
group_by(site) %>%
summarise(shan=diversity(tot_count)) %>%
ungroup() %>%
mutate(shan=ifelse(shan==0, min(shan[shan!=0]), shan))
# # calc shannon waterbirds
# shan_div_wb <- abundance %>% left_join(attributes) %>% left_join(inddf) %>%
# filter(Group=="waterbird") %>%
# # select(farm, pc_station_id, month_id, func_spec, max_count) %>%
# group_by(site, func_spec) %>%
# summarise(tot_count=sum(mean)) %>%
# ungroup() %>%
# group_by(site) %>%
# summarise(shan_wb=diversity(tot_count)) %>%
# ungroup() %>%
# mutate(shan_wb=ifelse(shan_wb==0, min(shan_wb[shan_wb!=0]), shan_wb))
# DT::datatable(shan_div)Once the components are done, BFI are calculated for each site and presented as scaled by plogis and scale with mean = 0
# merge to bfi table
bfi_site <- score_wt_counts %>% left_join(shan_div) %>%
mutate(bfi_unscaled=tot_wt_count*shan, bfi_scaled_plogis=plogis(as.numeric(scale(bfi_unscaled))),
bfi_scaled_scale=scale(as.numeric(scale(bfi_unscaled))) ) # %>%
#left_join(score_wt_counts_wb) %>%
#left_join(shan_div_wb) %>%
# rename(month=month_id, station=pc_station_id) %>%
# mutate(bfi_wb_unscaled=tot_wt_count_wb*shan_wb,
# bfi_wb_scaled=plogis(as.numeric(scale(bfi_wb_unscaled))))
# see table
# kbl(bfi_site)
# kable(bfi_site)
# see table
gt(bfi_site) %>%
tab_header(
title = md("Bird Friendliness Index **Parita sites**")#,
# subtitle = md("`gtcars` is an R dataset")
) # |> opt_interactive() | Bird Friendliness Index Parita sites | |||||
|---|---|---|---|---|---|
| site | tot_wt_count | shan | bfi_unscaled | bfi_scaled_plogis | bfi_scaled_scale |
| Parita_01_BayanoAvispa_03 | 1332.923 | 0.5495110 | 732.4557 | 0.21873900 | -1.27303004 |
| Parita_01_BayanoLasTablas_04 | 1803.923 | 0.5250362 | 947.1248 | 0.77322659 | 1.22662081 |
| Parita_01_ElAgallito_05 | 1715.105 | 0.5301453 | 909.2549 | 0.68689790 | 0.78565633 |
| Parita_01_ElReten_07 | 1657.773 | 0.5158807 | 855.2128 | 0.53901518 | 0.15637862 |
| Parita_01_EsteroCamaron_13 | 1587.714 | 0.5497024 | 872.7705 | 0.58923974 | 0.36082339 |
| Parita_01_HondaSalada_01 | 1845.604 | 0.5071854 | 936.0635 | 0.74985155 | 1.09782073 |
| Parita_01_LasPalmitas_14 | 1411.233 | 0.6011166 | 848.3156 | 0.51900731 | 0.07606591 |
| Parita_01_LosAzules_12 | 1575.255 | 0.5304803 | 835.6415 | 0.48212925 | -0.07151345 |
| Parita_01_LosPalacios_06 | 1693.108 | 0.5341873 | 904.4367 | 0.67470689 | 0.72955171 |
| Parita_01_LosPalacios_10 | 1823.275 | 0.5029117 | 916.9465 | 0.70583043 | 0.87521868 |
| Parita_01_LosPalacios_11 | 1716.573 | 0.5219536 | 895.9714 | 0.65271166 | 0.63098007 |
| Parita_01_Mamoni_08 | 1612.381 | 0.5114779 | 824.6970 | 0.45042486 | -0.19895422 |
| Parita_01_Nanzal_09 | 1265.006 | 0.5515041 | 697.6560 | 0.15732795 | -1.67824534 |
| Parita_01_PenonHonda_02 | 1548.695 | 0.5405983 | 837.2217 | 0.48672468 | -0.05311374 |
| Parita_01_Sarigua_15 | 1510.339 | 0.5352787 | 808.4522 | 0.40417207 | -0.38811122 |
| Parita_01_Sarigua_16 | 1146.096 | 0.5639216 | 646.3083 | 0.09311771 | -2.27614824 |