Abundance estimates for all species in Parita Bay, Panama
Site abundance prediction per species
Authors
Diego J. Lizcano
Jorge Velásquez-Tibata
Published
May 17, 2024
Fit the occupancy model of Royle and Nichols (2003), which relates probability of detection of the species to the number of individuals available for detection at each site.
Royle, J. A. and Nichols, J. D. (2003) Estimating Abundance from Repeated Presence-Absence Data or Point Counts. Ecology, 84(3) pp. 777–790.
The number of animals available for detection at site i is modelled as Poisson: \[
N_1 \sim Poisson(\lambda_i)
\]
We assume that all individuals (of same species) at site i during sample j have identical detection probabilities, \(r_{i,j}\) and and that detections are independent. The species will be recorded if at least one individual is detected. Thus, the detection probability for the species is linked to the detection probability for an individual by:
\[
p_{i,j} = 1-(1-r_{i,j})^{N_i}
\] The equation for the detection history is then: \[
y_{i,j} \sim Bernoulli(p_{i,j})
\] We used the unmarked package to make the estimates.
Read data
We read the data from excel and csv.
We prepared the data making a list of detection histories per species using the custom function f.matrix.creator2 and later collapsing the detection histories to 3 days using the custom function f.shrink.matrix.h.to9
Code
# load the custom functions to create a list of detection histories per species source("C:/CodigoR/AudubonPanama/R/matrix_creator.R")############### load Packageslibrary(readr) # readr library(readxl) # read excellibrary(plyr) library(dplyr)library(unmarked) # occupancy by maximum likelihoodlibrary(sf) # simple feature to get coordslibrary(DT) # to make data table with buttonslibrary(terra) # to convert sf to vectorlibrary(tidyverse) # some times ############### load data# old # parita_data_full <- read_delim("C:/CodigoR/AudubonPanama/data/all_records_audubon_panama_to_occu.csv", # delim = "\t", escape_double = FALSE, # col_types = cols(eventDate = col_date(format = "%Y-%m-%d")), # trim_ws = TRUE)# newparita_data_full <-read_excel("C:/CodigoR/AudubonPanama/data/all_records_audubon_panama_Darwin_Core_revision.xlsx", col_types =c("text", "text", "text", "date", "text", "text", "numeric", "numeric", "text", "text", "numeric", "text", "numeric", "numeric", "text", "numeric", "numeric", "numeric", "numeric", "numeric", "text", "text", "numeric", "text", "numeric", "text", "text", "text", "text", "text", "numeric", "text", "numeric", "text", "numeric", "text", "text"))covs <-read_csv("C:/CodigoR/AudubonPanama/shp/sites_covs_parita_nona.csv")## Nyctidromus albicollis## Aramus guarauna## Buteogallus anthracinus## Dryocopus lineatus## Piranga rubra# filter by threshold and speciesparita_data_allsp <- parita_data_full %>%filter(confidence >= threshold)# %>%# filter(scientificName == c("Nyctidromus albicollis", # "Aramus guarauna",# "Buteogallus anthracinus",# "Dryocopus lineatus",# "Piranga rubra"# )) # insert Max and min to get matrix limits# start_dateparita_data2 <- parita_data_allsp %>%group_by(locationID) %>%mutate(start_date =min(eventDate),end_date=max(eventDate))# apply function to get matrix per speciesparita <-f.matrix.creator2(data = parita_data2, year =2023)# # centroide for terra# # centoide <- centroids(puntos, TRUE)# centroide <- c(mean(as.numeric(parita_data_allsp$decimalLongitude)), mean(as.numeric(parita_data_allsp$decimalLatitude)))# clip_window <- extent(-80.67105, -80.05822, 7.713744, 8.489888) # 7.932311, -80.612233 8.360073, -80.180109# bb <- c(-80.612233, -80.180109, 7.932311, 8.360073)# srtm <- raster::getData('SRTM', lon=centroide[1], lat=centroide[2])# crop the raster using the vector extent# srtm_crop <- raster::crop(srtm, clip_window)#rast(srtm_crop) # convert to terra# altitud <- elevation_3s(-72.893262, 7.664081007, path="data")# altitud <- as.numeric(Camptostoma_obsoletum_occu_dia$Altitude)# convierte covs a puntos terrapuntos <-vect(covs, geom=c("Longitude", "Latitude"), crs="EPSG:4326")# convierte a sfpuntos_sf <- sf::st_as_sf(puntos)######### precipitacionprecip_Anton <-read_csv("C:/CodigoR/AudubonPanama/data/Precip_Anton.csv")precip <-data.frame(f.shrink.matrix.h.to9(precip_Anton))precip[,1] <-0#fix start#
Function to get abundance
These methods return an object storing the posterior distributions of the latent variables at each site. We used the empirical Bayes methods used by unmarked::ranef which can underestimate the variance of the posterior distribution because they do not account for uncertainty in the hyperparameters (lambda). Note also that the posterior mode appears to exhibit some bias as an estimator or abundance. Consider using the posterior mean instead, even though it will not be an integer in general.
One species, single season model in a custom function
The function ranef_by_sp calculates the abundance per species.
Code
ranef_by_sp <-function(spx=1){ suppressWarnings({ #species name# Make UMF objectumf_y_full_spx<-unmarkedFrameOccu(y=f.shrink.matrix.to9( # y collapsed to 3daysas.data.frame(parita$occur[[spx]])),siteCovs=NULL,obsCovs=list(prec=precip,hour=data.frame(hour =f.shrink.matrix.h.to9(parita$dete[[spx]])) ) ) # fit unmarkedfit_unm_spx_0 <- unmarked::occuRN(~1~1, data=umf_y_full_spx) # ok!fit_unm_spx_p <- unmarked::occuRN(~scale(prec) ~1, data=umf_y_full_spx) # ok!fit_unm_spx_h <- unmarked::occuRN(~scale(hour) ~1, data=umf_y_full_spx) # ok!# fit list for occupancyfms2<-fitList("p(.) lambda(.)"=fit_unm_spx_0,"p(prec) lambda(.)"=fit_unm_spx_p,"p(hour) lambda(.)"=fit_unm_spx_h)# #"p(h) Ocu(.)"=fit_unm_sp2_0h, #,# #"p(c) Ocu(.)"=fit_unm_sp2_0c,# "p(.) lambda(BGB_Spawn)"=fit_unm_sp2_1,# "p(.) lambda(human_foot)"=fit_unm_sp2_2,# "p(.) lambda(NDVI)"=fit_unm_sp2_3,# "p(.) lambda(river)"=fit_unm_sp2_4,# #"p(.) Ocu(river^2)"=fit_unm_sp2_4s,# "p() lambda(carbon)"=fit_unm_sp2_5,# "p() lambda(canopy)"=fit_unm_sp2_6,# "p() lambda(canopy+BGB_Spawn)"=fit_unm_sp2_7,# "p() lambda(BGB_Spawn+human_foot)"=fit_unm_sp2_8,# "p() lambda(ALL)"=fit_unm_sp2_9# )# model selection detection unmarked models <-modSel(fms2) models# Empirical Bayes estimates of abundance at each site by speciesre_0 <-ranef(fit_unm_spx_0)re_h <-ranef(fit_unm_spx_h)# plot(re)# post.df <- as(re_0, "data.frame")###################################################### to save the ranef result (S3 Object) as dataframe ######################################################abund_per_site <-data.frame(sp=names(parita$dete[spx]), site=names( parita$occu[[1]][,1]),mode=(bup(re_0, stat="mode") ),mean=(bup(re_0, stat="mean") ), confint(re_0, level=0.9))}) # no warning part# return(list(paste("###", sp=names(parita$dete[spx]) ), # sp name# models, # model AICs# Null=re_0, # mull model estimates# hour=re_h)) # hour in detection estimates# return to export tablereturn(abund_per_site)} # end of function
use the function
Runs a loop on a list containing the species, and later making a table of species sites and abundances.
Table of abundances
Code
# library (plyr)# Convert the list to a dataframedf <-ldply(abun_site, data.frame)# save to table# write.csv(df, "C:/CodigoR/AudubonPanama/data/abundance_result.csv")# make table with buttons datatable( df, extensions ='Buttons', options =list(dom ='Bfrtip',buttons =c('copy', 'csv') ))
