Code
library(leaflet)
library(ggplot2)
library(vegan)
library(gt)
library(knitr)
library(kableExtra)
library(envalysis)
library(sads)
library(effects)
library(unmarked)Bird species richness, abundance and density using point count method
Abstract
A short tutorial to explore bird species richness, distribution and abundance
Keywords
bird, point count method, species richness, abundance
Bird point count method
Summary (according to google AI): The bird point count method is a widely used, standardized technique where an observer stands at a fixed location (a “point”) for a set time (e.g., 3-10 mins), recording all birds detected by sight or sound to estimate species richness and relative abundance, crucial for monitoring bird populations and habitats. Observers typically count birds within a certain radius (or sometimes unlimited distance), noting species, count, and sometimes distance, often early in the morning during breeding seasons, using tools like binoculars and datasheets, and recording environmental conditions like wind (Fontúrbel 2020).
Key Steps
Fixed Location: An observer stays at a single spot (station) for a predetermined duration.
Standardized Time: Counts usually last 3 to 10 minutes, with variations for different habitats.
Detection: All birds seen or heard (singing, calls) are tallied.
Distance: Some protocols use a fixed radius (e.g., 100m) or record sighting distances to estimate density; others use unlimited distance.
Environmental Data: Wind speed, cloud cover, and habitat type are recorded.
Timing: Often done in the morning during peak singing (early breeding season) under good weather.
Purpose: To assess bird distribution, abundance, and habitat health over time.
Methods
Select Points: Choose stations (often GPS located) spread out enough to avoid double-counting birds.
Record Data: At each point, record species, numbers, and potentially distance for a set time.
Analyze: Data is aggregated to understand trends in bird populations, helping managers assess conservation success.
eBird helps to create checklist and save your record species: How to count birds (eBird)
Read the datasets
Download count data and spatial data sets
Load R libraries
Load datasets
Location and detection distance
Map the location of the point counts
buffer (100 meters) and species spotted (checklist in eBird)
You can view here one of the station uploaded in the eBird checklist
Code
leaflet() %>%
addTiles() %>%
addCircleMarkers(data = bird_data,
label = ~as.character(Species),
popup = leafpop::popupTable(bird_data),
fillColor = "blue", weight = 0.5,
~ long,
~ lat,
group = "birds",
clusterOptions = markerClusterOptions(),
stroke = FALSE, fillOpacity = 0.6) %>%
addPolygons(data= poly_sampling, color = "orange", weight = 3, smoothFactor = 0.5,
opacity = 1.0, fillOpacity = 0.01,
group="sampling area") %>%
addPolygons(data= buffer_sampling, color = "blue", weight = 1, smoothFactor = 0.5,
opacity = 1.0, fillOpacity = 0.01,
group="buffer") %>%
addProviderTiles(providers$Esri.WorldStreetMap,group ="WSM") %>%
addTiles(group = "OSM") %>%
addTiles(urlTemplate = "https://mts1.google.com/vt/lyrs=s&hl=en&src=app&x={x}&y={y}&z={z}&s=G",
attribution = 'Google',group ="Google") %>%
addLayersControl(baseGroups = c("WSM","OSM","Google"),
overlayGroups = c("birds","sampling area","buffer"),
options = layersControlOptions(collapsed = FALSE)) %>%
addScaleBar(position = "bottomright",
scaleBarOptions(maxWidth = 100, metric = TRUE, imperial = FALSE,
updateWhenIdle = TRUE)) Assess the distance (100 meters)
Exemple using the records of the Asian Openbill (Anastomus oscitans)
Code
#|
bird_count_data_select <- distance_data |>
dplyr::filter(species == "Asian Openbill (Anastomus oscitans)")
ggplot(bird_count_data_select, aes(x =distance)) +
geom_density() +
theme_minimal() +
labs(title = "Density plot of Asian Openbill in relation to distance",
y = "Density",
x = "Distance (in meters)",
color="orange") +
theme_publish()
Exemple using the records of the Baya Weaver (Ploceus philippinus)
Code
#|
distance_data |>
dplyr::filter(species == "Baya Weaver (Ploceus philippinus)") |>
ggplot(aes(x =distance)) +
geom_density() +
theme_minimal() +
labs(title = "Density plot of Baya Weaver in relation to distance",
y = "Density",
x = "Distance (in meters)",
color="orange") +
theme_publish()
Species richness
Species richness by location
Code
#|
bird_data |>
dplyr::group_by(eBird_checklist, Species) |>
dplyr::summarise(total = sum(Count), .groups = "drop") |>
dplyr::group_by(eBird_checklist) |>
dplyr::summarise(richness = dplyr::n_distinct(Species)) |>
ggplot(aes(x = eBird_checklist, y = richness, fill = eBird_checklist)) +
geom_col() +
labs(y = "Observed Species Richness", x = "Location") +
theme_publish() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
theme(legend.position = "none")
Code
#|
bird_data |>
dplyr::group_by(eBird_checklist, Species) |>
dplyr::summarise(total = sum(Count), .groups = "drop") |>
dplyr::group_by(eBird_checklist) |>
dplyr::summarise(richness = dplyr::n_distinct(Species)) |>
dplyr::left_join(bird_data |>
dplyr::select(eBird_checklist, long,lat),
by = "eBird_checklist") |>
leaflet() %>%
addTiles() %>%
addCircleMarkers(
lng = ~long,
lat = ~lat,
radius = ~ richness, # Scale the values (e.g., using sqrt and division for better visual scaling)
popup = ~as.character(richness),
stroke = FALSE,
fillOpacity = 0.5,
group="species richness") %>%
addPolygons(data= poly_sampling, color = "orange", weight = 3, smoothFactor = 0.5,
opacity = 1.0, fillOpacity = 0.01,
group="sampling area") %>%
addPolygons(data= buffer_sampling, color = "red", weight = 1, smoothFactor = 0.5,
opacity = 1.0, fillOpacity = 0.01,
group="buffer") %>%
addProviderTiles(providers$Esri.WorldStreetMap,group ="WSM") %>%
addTiles(group = "OSM") %>%
addTiles(urlTemplate = "https://mts1.google.com/vt/lyrs=s&hl=en&src=app&x={x}&y={y}&z={z}&s=G",
attribution = 'Google',group ="Google") %>%
addLayersControl(baseGroups = c("WSM","OSM","Google"),
overlayGroups = c("species richness","sampling area","buffer"),
options = layersControlOptions(collapsed = FALSE)) %>%
addScaleBar(position = "bottomright",
scaleBarOptions(maxWidth = 100, metric = TRUE, imperial = FALSE,
updateWhenIdle = TRUE))Species richness estimators
Create a species-by-site matrix (rows = sites, cols = species)
Code
#|
comm_matrix <- bird_data |>
dplyr::group_by(eBird_checklist, Species) |>
dplyr::summarise(Count = sum(Count), .groups = "drop") |>
tidyr::pivot_wider(names_from = Species, values_from = Count, values_fill = 0) |>
tibble::column_to_rownames("eBird_checklist")
comm_matrix |>
gt(
)| Asian Openbill (Anastomus oscitans) | Black Drongo (Dicrurus macrocercus) | Black-naped Oriole (Oriolus chinensis) | Coppersmith Barbet (Psilopogon haemacephalus) | Great Egret (Ardea alba) | Great Myna (Acridotheres grandis) | House Swift (Apus nipalensis) | Indochinese Roller (Coracias affinis) | Lesser Whistling-Duck (Dendrocygna javanica) | Oriental Darter (Anhinga melanogaster) | Red-wattled Lapwing (Vanellus indicus) | Rock Pigeon (Feral Pigeon) (Columba livia (Feral Pigeon)) | Scarlet-backed Flowerpecker (Dicaeum cruentatum) | Siamese Pied Starling (Gracupica floweri) | Common Myna (Acridotheres tristis) | Javan Pond-Heron (Ardeola speciosa) | Malayan Night Heron (Gorsachius melanolophus) | Ornate Sunbird (Cinnyris ornatus) | Purple Heron (Ardea purpurea) | Streak-eared Bulbul (Pycnonotus conradi) | Blue-tailed Bee-eater (Merops philippinus) | Indian Cormorant (Phalacrocorax fuscicollis) | Little Grebe (Tachybaptus ruficollis) | Scaly-breasted Munia (Lonchura punctulata) | Spotted Owlet (Athene brama) | White-breasted Waterhen (Amaurornis phoenicurus) | Asian Koel (Eudynamys scolopaceus) | Glossy Ibis (Plegadis falcinellus) | Large-billed Crow (Corvus macrorhynchos) | Lineated Barbet (Psilopogon lineatus) | Oriental Magpie-Robin (Copsychus saularis) | Painted Stork (Mycteria leucocephala) | Pied Kingfisher (Ceryle rudis) | Baya Weaver (Ploceus philippinus) | Chinese Pond-Heron (Ardeola bacchus) | Plain Prinia (Prinia inornata) | Zebra Dove (Geopelia striata) | Asian Brown Flycatcher (Muscicapa dauurica) | Gray-headed Swamphen (Porphyrio poliocephalus) | White-throated Kingfisher (Halcyon smyrnensis) | Black-browed Reed Warbler (Acrocephalus bistrigiceps) | Brown Shrike (Lanius cristatus) | Chestnut-headed Bee-eater (Merops leschenaulti) | Common Kingfisher (Alcedo atthis) | Common Sandpiper (Actitis hypoleucos) | Paddyfield Pipit (Anthus rufulus) | Yellow-bellied Prinia (Prinia flaviventris) | Barn Swallow (Hirundo rustica) | Little Cormorant (Microcarbo niger) | Spotted Dove (Spilopelia chinensis) | Cotton Pygmy-Goose (Nettapus coromandelianus) | Black Baza (Aviceda leuphotes) | Yellow Bittern (Botaurus sinensis) | Yellow-vented Bulbul (Pycnonotus goiavier) | White-browed Crake (Poliolimnas cinereus) | Green-billed Malkoha (Phaenicophaeus tristis) | Little Egret (Egretta garzetta) | Medium Egret (Ardea intermedia) | Taiga Flycatcher (Ficedula albicilla) | Black-winged Stilt (Himantopus himantopus) | Eastern Cattle-Egret (Ardea coromanda) | Gray-backed Shrike (Lanius tephronotus) | Ashy Woodswallow (Artamus fuscus) | Shikra (Tachyspiza badia) | Red Collared-Dove (Streptopelia tranquebarica) | Gray Heron (Ardea cinerea) |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | 1 | 2 | 1 | 1 | 23 | 13 | 1 | 23 | 2 | 2 | 2 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 14 | 2 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 1 | 0 | 3 | 1 | 0 | 2 | 4 | 1 | 1 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 13 | 0 | 0 | 0 | 3 | 0 | 3 | 0 | 0 | 0 | 3 | 5 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 3 | 2 | 1 | 2 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 34 | 0 | 0 | 0 | 0 | 6 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 2 | 0 | 0 | 1 | 0 | 2 | 0 | 0 | 0 | 1 | 0 | 1 | 1 | 27 | 1 | 1 | 1 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 10 | 2 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 26 | 0 | 0 | 0 | 0 | 0 | 0 | 2 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 112 | 1 | 4 | 3 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 7 | 0 | 1 | 0 | 2 | 0 | 0 | 0 | 9 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 2 | 2 | 0 | 0 | 0 | 2 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 200 | 0 | 2 | 0 | 1 | 2 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 10 | 1 | 0 | 0 | 0 | 0 | 5 | 0 | 11 | 0 | 2 | 1 | 0 | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 2 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 5 | 0 | 3 | 0 | 0 | 2 | 0 | 2 | 1 | 5 | 1 | 1 | 1 | 2 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 3 | 2 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 11 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 2 | 0 | 0 | 0 | 1 | 3 | 1 | 0 | 0 | 0 | 0 | 2 | 1 | 2 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 2 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 4 | 1 | 0 | 4 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 32 | 0 | 2 | 1 | 0 | 0 | 1 | 1 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 3 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 5 | 2 | 0 | 0 | 1 | 0 | 0 | 0 | 10 | 0 | 0 | 3 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 150 | 0 | 0 | 0 | 0 | 0 | 0 | 3 | 1 | 1 | 0 | 0 | 0 | 0 | 4 | 1 | 0 | 0 | 4 | 1 | 3 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 16 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 4 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 11 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 2 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 1 | 2 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 2 | 1 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 3 | 1 | 0 | 1 | 2 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 1 | 0 | 1 | 2 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 1 | 2 | 0 | 0 | 6 | 7 | 0 | 1 | 0 | 0 | 0 | 2 | 0 | 1 | 0 | 5 | 0 | 0 | 0 | 2 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 4 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 10 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 10 | 6 | 0 | 7 | 35 | 1 | 0 | 0 | 0 | 0 |
| 8 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 3 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 2 | 0 | 4 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 1 | 1 | 0 | 0 |
| 2 | 1 | 0 | 0 | 0 | 0 | 3 | 0 | 0 | 0 | 0 | 2 | 0 | 0 | 2 | 0 | 0 | 0 | 0 | 0 | 2 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 2 | 0 | 0 | 0 | 0 | 2 | 2 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 1 | 0 |
| 19 | 9 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 3 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 39 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 3 | 0 | 0 | 0 |
| 5 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 3 | 2 | 2 | 0 | 0 | 0 | 0 | 1 | 1 | 0 | 2 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 3 | 0 | 0 | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 12 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 2 | 0 | 0 | 1 |
Estimate richness using Chao1 and other estimators
Code
#|
pool <- specpool(comm_matrix)
kable(pool, caption = "Species Pool Estimators") %>%
kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))| Species | chao | chao.se | jack1 | jack1.se | jack2 | boot | boot.se | n | |
|---|---|---|---|---|---|---|---|---|---|
| All | 66 | 99.88235 | 19.20208 | 88.58824 | 6.331117 | 103.1471 | 75.89184 | 3.472629 | 17 |
Code
spec_accum <- specaccum(comm_matrix, method = "random")
plot(spec_accum, ci.type = "poly", col = "blue", lwd = 2,
ci.lty = 0, ci.col = "lightblue",
xlab = "Number of Locations", ylab = "Species Richness")
Abundance
Abundance of bird species per location
Abundance for species with more than 10 individuals counted
Code
#|
ggplot(bird_data|>
dplyr::filter(Count > 10),
aes(x = eBird_checklist, y = Count, fill = Species)) +
geom_col() +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
labs(
title = "Bird abundance by location",
x = "Location",
y = "Count"
)
Abundance distribution models
Abundance distribution models (ADMs), and Species Abundance Distributions (SADs), describe how many individuals of each species exist in a community, revealing patterns of commonness (few species) and rarity (many species), often forming a “hollow curve”. Explaining them involves showing this characteristic shape and linking it to underlying ecological theories like niche partitioning (different resources) or neutral processes (random chance), with newer approaches like Powerbend.
Key models include the logseries (many rare), lognormal (internal mode), and broken-stick (niche-based), helping understand biodiversity, community assembly, and responses to disturbances. (Gao 2025)
The “Hollow Curve”: A graph showing the number (or % of species) on the y-axis and abundance (number of individuals) on the x-axis, which almost always shows a few very abundant species and many rare ones, forming a J-shape or hollow curve, which seems universal in nature.
Some models Logseries (Fisher): Predicts many singletons (species with only one individual) and no internal peak.
Lognormal (Preston): More realistic; reveals a peak of intermediate abundance when rare species are accounted for, creating a bell-curve when plotted on a log scale.
Broken-Stick (MacArthur): a Niche-Based model (ecological mechanisms), with the species divide available resources (niches) randomly, like breaking a stick.
Neutral Models: Assume all species are ecologically identical, with differences driven by random birth, death, and dispersal.
Pueyo’s Power-bend model: a discrete probability distribution two-parameter model that generalizes the power-law and log-series distributions; specifically designed to correct the tendency of power-law models to overestimate the abundance of very common species in a sample.
Distribution of bird species abundance (“Hollow Curve”)
Code
#|
# 1 Aggregate counts by species (sum across locations)
sp_abund <- bird_data |>
dplyr::group_by(Species) |>
dplyr::summarise(TotalCount = sum(Count, na.rm = TRUE)) |>
dplyr::arrange(desc(TotalCount))
# 2 Convert to numeric abundance vector
abund_vector <- sp_abund$TotalCount
# 3 Plot species rank-abundance curve
rank_data <- data.frame(
Rank = 1:length(abund_vector),
Abundance = sort(abund_vector, decreasing = TRUE)
)
ggplot(rank_data, aes(x = Rank, y = Abundance)) +
geom_point(size = 2) +
geom_line() +
scale_y_log10() +
theme_minimal(base_size = 14) +
labs(title = "Species Rank–Abundance Curve",
x = "Species Rank",
y = "Abundance (log scale)") +
theme_publish() 
Fit distribution models
Code
#|
# Fit logseries
fit_logseries <- fitsad(abund_vector, "ls")
# Fit lognormal
fit_lnorm <- fitsad(abund_vector, "lnorm")
# Fit broken-stick
fit_bs <- fitsad(abund_vector, "bs")
# Fit a neutral model
fit_neutral <- fitsad(abund_vector, sad = "volkov")
# Fit Pueyo's Power-bend model
fit_powerbend <- fitpowbend(abund_vector)
# Compute AICs with names
aic_values <- AIC(
fit_logseries,
fit_lnorm,
fit_bs,
fit_neutral,
fit_powerbend
)
labels <- c(
fit_logseries = "Log-series",
fit_lnorm = "Lognormal",
fit_bs = "Broken-stick",
fit_neutral = "Neutral theory",
fit_powerbend = "Pueyo's Power-bend"
)
# Convert to a data frame
aic_df <- data.frame(
Model = labels,
AIC = aic_values$AIC
)
# Rank models (1 = best)
aic_df$Rank <- rank(aic_df$AIC, ties.method = "first")
# Sort by rank
aic_df <- aic_df[order(aic_df$Rank), ]
aic_df |>
gt(
)| Model | AIC | Rank |
|---|---|---|
| Pueyo's Power-bend | 446.8116 | 1 |
| Log-series | 447.3464 | 2 |
| Neutral theory | 449.3651 | 3 |
| Lognormal | 456.9042 | 4 |
| Broken-stick | 535.7691 | 5 |
Code
# Extract coefficients table
model_pwb <- summary(fit_powerbend)
coef_tbl <- as.data.frame(model_pwb@coef)
coef_tbl$term <- rownames(coef_tbl)
coef_tbl <- coef_tbl %>% dplyr::relocate(term)
knitr::kable(coef_tbl, caption = "Pueyo's Power-bend Model Coefficients")| term | Estimate | Std. Error | z value | Pr(z) | |
|---|---|---|---|---|---|
| s | s | 1.186476 | 0.1138346 | 10.42281 | 0 |
| oM | oM | 5.036158 | 0.5847058 | 8.61315 | 0 |
Code
# statistics
stats_tbl <- tibble::tibble(
logLik = as.numeric(logLik(fit_powerbend)),
AIC = AIC(fit_powerbend),
BIC = BIC(fit_powerbend),
n = length(fit_powerbend@data),
model = "Power-bend"
)
knitr::kable(stats_tbl, caption = "Pueyo's Power-bend Model Fit Statistics")| logLik | AIC | BIC | n | model |
|---|---|---|---|---|
| -221.4058 | 446.8116 | 451.1909 | 1 | Power-bend |
Plot the best model (here Pueyo’s Power-bend)
Code
#|
# Disable "Press Return to see next plot"
op <- par(ask = FALSE)
par(mfrow = c(2, 2))
plot(fit_powerbend, main = "Pueyo's Power-bend fit to SAD",
col = "darkblue", lwd = 2)
Individual abundance and density
First method using general linear modeling (GLM)
Exemple using the records of the Asian Openbill (Anastomus oscitans)
Code
#|
data_AsianOpenbill <- bird_data |>
dplyr::filter(Species == "Asian Openbill (Anastomus oscitans)")
m1 <- glm(Count ~ long + lat,
data = data_AsianOpenbill,
family = poisson)
# Coefficients
tidy_m1 <- broom::tidy(m1, conf.int = TRUE)
# Model statistics
stats_m1 <- broom::glance(m1) |>
dplyr::select(AIC, deviance, df.residual, logLik)
# TABLE 1: coefficients
gt(tidy_m1) |>
tab_header(title = "Poisson GLM Coefficients") |>
fmt_number(columns = where(is.numeric), decimals = 3)| Poisson GLM Coefficients | ||||||
|---|---|---|---|---|---|---|
| term | estimate | std.error | statistic | p.value | conf.low | conf.high |
| (Intercept) | 1,219.071 | 2,280.516 | 0.535 | 0.593 | −3,300.824 | 5,648.559 |
| long | −21.461 | 23.002 | −0.933 | 0.351 | −66.323 | 23.953 |
| lat | 59.620 | 22.537 | 2.645 | 0.008 | 15.739 | 104.135 |
Code
# TABLE 2: model fit
gt(stats_m1) |>
tab_header(title = "Model Fit Statistics") |>
fmt_number(columns = where(is.numeric), decimals = 3)| Model Fit Statistics | |||
|---|---|---|---|
| AIC | deviance | df.residual | logLik |
| 180.781 | 114.468 | 14.000 | −87.390 |
Code
effects = allEffects(m1)
plot(effects,
col = 2,
ylab = "Number (residuals)",
type = "response")
Estimated abundance by location
Code
#|
data_AsianOpenbill$predicted_count <- predict(
m1,
type = "response" # gives predicted counts, NOT log-scale
)
ggplot(data_AsianOpenbill, aes(x = pointname, y = predicted_count)) +
geom_col() +
labs(
title = "Predicted Count per Point",
x = "Point",
y = "Predicted Count"
) +
theme_publish() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) 
Map
Code
#|
leaflet(data = data_AsianOpenbill) %>%
addTiles() %>%
addCircleMarkers(
lng = ~long,
lat = ~lat,
radius = ~ Count,
# popup = leafpop::popupTable(bird_data),
stroke = FALSE,
fillOpacity = 0.5,
group="observed abundance") %>%
addCircleMarkers(
lng = ~long,
lat = ~lat,
radius = ~ as.integer(round(predicted_count)),
color = "red",
popup = leafpop::popupTable(data_AsianOpenbill),
stroke = FALSE,
fillOpacity = 0.5,
group="predicted abundance") %>%
addPolygons(data= poly_sampling, color = "orange", weight = 3, smoothFactor = 0.5,
opacity = 1.0, fillOpacity = 0.01,
group="sampling area") %>%
addPolygons(data= buffer_sampling, color = "red", weight = 1, smoothFactor = 0.5,
opacity = 1.0, fillOpacity = 0.01,
group="buffer") %>%
addProviderTiles(providers$Esri.WorldStreetMap,group ="WSM") %>%
addTiles(group = "OSM") %>%
addTiles(urlTemplate = "https://mts1.google.com/vt/lyrs=s&hl=en&src=app&x={x}&y={y}&z={z}&s=G",
attribution = 'Google',group ="Google") %>%
addLayersControl(baseGroups = c("WSM","OSM","Google"),
overlayGroups = c("observed abundance","predicted abundance","sampling area","buffer"),
options = layersControlOptions(collapsed = FALSE)) %>%
addScaleBar(position = "bottomright",
scaleBarOptions(maxWidth = 100, metric = TRUE, imperial = FALSE,
updateWhenIdle = TRUE)) Density and total abundance
First estimation of density and total abundance
Exemple using the records of the Asian Openbill (Anastomus oscitans)
Code
#|
data_AsianOpenbill$density = data_AsianOpenbill$predicted_count / (pi * 100 ^2)
abundance_mean <- mean(data_AsianOpenbill$predicted_count)
density_mean <- mean(data_AsianOpenbill$density)
total_area = 2326317
total_number <- mean(data_AsianOpenbill$density) * total_area
cbind.data.frame(abundance_mean, density_mean,total_number) |>
dplyr::rename(`abundance (mean among plot)`= abundance_mean,
`density (square meters)` = density_mean,
`total estimated abundance (total sampling area)` = total_number) |>
gt()| abundance (mean among plot) | density (square meters) | total estimated abundance (total sampling area) |
|---|---|---|
| 8.823529 | 0.0002808617 | 653.3733 |
Second estimation of density (using unmarked)
Exemple using the records of the Asian Openbill (Anastomus oscitans)
Code
# Pivot to wide format
birds_wide <- bird_data |>
dplyr::filter(Species == "Asian Openbill (Anastomus oscitans)") |>
dplyr::select(eBird_checklist, pointname, Count) |>
tidyr::pivot_wider(names_from = pointname, values_from = Count,
names_prefix = "count_visit_")
birds_wide |>
gt()| eBird_checklist | count_visit_P001 | count_visit_P002 | count_visit_P003 | count_visit_P004 | count_visit_P005 | count_visit_P006 | count_visit_P007 | count_visit_P008 | count_visit_P009 | count_visit_P010 | count_visit_P011 | count_visit_P012 | count_visit_P013 | count_visit_P014 | count_visit_P015 | count_visit_P016 | count_visit_P017 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| S282429469 | 1 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| S282434325 | NA | 14 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| S282436054 | NA | NA | 13 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| S282438024 | NA | NA | NA | 34 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| S282439163 | NA | NA | NA | NA | 10 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| S282440354 | NA | NA | NA | NA | NA | 7 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| S282441414 | NA | NA | NA | NA | NA | NA | 10 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| S282452204 | NA | NA | NA | NA | NA | NA | NA | 3 | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| S282454221 | NA | NA | NA | NA | NA | NA | NA | NA | 1 | NA | NA | NA | NA | NA | NA | NA | NA |
| S282455845 | NA | NA | NA | NA | NA | NA | NA | NA | NA | 5 | NA | NA | NA | NA | NA | NA | NA |
| S282458258 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | 16 | NA | NA | NA | NA | NA | NA |
| S282462040 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | 1 | NA | NA | NA | NA | NA |
| S282466316 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | 1 | NA | NA | NA | NA |
| S282468198 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | 8 | NA | NA | NA |
| S282469766 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | 2 | NA | NA |
| S282471178 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | 19 | NA |
| S282472982 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | 5 |
Code
# Convert to unmarkedFramePCount object
umf <- unmarkedFramePCount(
y = birds_wide |>
dplyr::select(starts_with("count_visit_")),
siteCovs = birds_wide |>
dplyr::select(eBird_checklist)
)
fm1 <- pcount(~1 ~1, data = umf, K = 50) # K = max expected count
# summary(fm1)
# Predict expected abundance per site
pred <- predict(fm1, type = "state")
birds_wide$predicted_abundance <- pred$Predicted
# Predict expected abundance per site
pred <- predict(fm1, type = "state")
birds_wide$predicted_abundance <- pred$PredictedCode
#|
data_AsianOpenbill$density2 = birds_wide$predicted_abundance/ (pi * 100 ^2)
predicted_abund <- mean(birds_wide$predicted_abundance)
density_mean <- mean(data_AsianOpenbill$density2)
total_area = 2326317
total_number <- mean(data_AsianOpenbill$density2) * total_area
cbind.data.frame(predicted_abund,density_mean,total_number) |>
dplyr::rename(`abundance (mean among plot)`= predicted_abund,
`density (square meters)` = density_mean,
`total estimated abundance (total sampling area)` = total_number) |>
gt()| abundance (mean among plot) | density (square meters) | total estimated abundance (total sampling area) |
|---|---|---|
| 14.39145 | 0.0004580942 | 1065.672 |
Resources
Using eBird Mobile to submit your sightings - eBird Essentials
References
Fontúrbel, Rodríguez-Gómez GB, FE. 2020. “Sampling Understory Birds in Different Habitat Types Using Point Counts and Camera Traps.” Ecological Indicators 119: 106863. https://doi.org/10.1016/j.ecolind.2020.106863.
Gao, Abdullah, Y. 2025. “The Powerbend Distribution Provides a Unified Model for the Species Abundance Distribution Across Animals, Plants and Microbes.” Nature Communications 16: 4035. https://doi.org/10.1038/s41467-025-59253-9.