This documents outlines how a single-taxon indicator (Aglais urticae on farms in Finland) can be calculated in R using the RBMS method.
The following packages are required. All packages are available on CRAN apart from
{fbi}
and {rbms}
which can be installed from
GitHub.
These five fields are required for the survey data.
These filters restrict the survey data to the “Butterflies in Finnish agricultural landscapes” dataset with the selected data fields have no missing data.
The survey data can now downloaded from FinBIF.
surveys <- finbif_occurrence(
filter = filter,
select = select,
aggregate = "events",
aggregate_counts = FALSE,
n = "all",
quiet = TRUE
)
A single processing function is applied to the survey data to convert the year month and day fields into a single date field.
Count data requires three fields to be selected: the survey
identifier (document_id
) the survey site section
(section
) and the measure of abundance
(abundance_interpreted
).
The count data requires the same filters as the survey data (though
the filter has_value
needs to be redefined).
The count data for Aglais urticae (Small tortoiseshell) can now be downloaded from FinBIF.
counts <- finbif_occurrence(
taxa = "Aglais urticae",
filter = filter,
select = select ,
n = "all",
quiet = TRUE
)
Two processing functions are applied to the count data to: sum the counts over the survey site sections; and combine the count and survey data together.
An RBMS model is used to estimate the change in abundance over time (see RBMS documentation for details) with uncertainty estimated via bootstrapping.
StartMonth <- 4
EndMonth <- 9
StartDay <- 1
EndDay <- 30
Anchor <- TRUE
AnchorLength <- 14
AnchorLag <- 14
AnchorTimeUnit <- "d"
NbrSample <- 500
MinVisit <- 5
MinOccur <- 2
MinNbrSite <- 5
MaxTrial <- 4
GamFamily <- "nb"
FlightCurveTimeUnit <- "w"
MultiVisit <- "mean"
MinFC <- 0.10
boot_n <- 200
surveys <- select(surveys, site_id = location_id, year, date)
init_year <- min(surveys[["year"]])
last_year <- max(surveys[["year"]])
nyears <- last_year - init_year + 1
ts_date <- ts_dwmy_table(InitYear = init_year, LastYear = last_year)
ts_season <- ts_monit_season(
d_series = ts_date,
StartMonth = StartMonth,
EndMonth = EndMonth,
StartDay = StartDay,
EndDay = EndDay,
Anchor = Anchor,
AnchorLength = AnchorLength,
AnchorLag = AnchorLag,
TimeUnit = AnchorTimeUnit
)
ts_season_visit <- ts_monit_site(m_visit = surveys, ts_season = ts_season)
counts <- mutate(counts, species = 1)
counts <- select(
counts, count = abundance, site_id = location_id, year, date, species
)
ts_season_count <- ts_monit_count_site(
m_season_visit = ts_season_visit, m_count = counts
)
ts_flight_curve <- flight_curve(
ts_season_count = ts_season_count,
NbrSample = NbrSample,
MinVisit = MinVisit,
MinOccur = MinOccur,
MinNbrSite = MinNbrSite,
MaxTrial = MaxTrial,
GamFamily = GamFamily,
SpeedGam = FALSE,
TimeUnit = FlightCurveTimeUnit,
MultiVisit = MultiVisit,
verbose = FALSE
)
impt_counts <- impute_count(
ts_season_count = ts_season_count,
ts_flight_curve = ts_flight_curve[["pheno"]],
TimeUnit = FlightCurveTimeUnit,
MultiVisit = MultiVisit
)
sindex <- site_index(butterfly_count = impt_counts, MinFC = MinFC)
bootsample <- boot_sample(data = sindex, boot_n = boot_n)
index_mc <- mapply(
collated_index,
bootID = seq_len(boot_n),
MoreArgs = list(data = sindex, s_sp = 1, boot_ind = bootsample),
SIMPLIFY = FALSE
)
An index of change in relative abundance with uncertainty derived from the bootstrapped estimates can be created by setting the base year to the year 2000.
base <- which(sort(unique(surveys[["year"]])) == 2000)
index_mc <- lapply(index_mc, getElement, "col_index")
index_mc <- do.call(rbind, index_mc)
index_mc <- mutate(index_mc, mc = log(pmax(1 / 100, COL_INDEX)), time = M_YEAR)
index_mc <- group_by(index_mc, BOOTi)
index_mc <- arrange(index_mc, time)
index_mc <- mutate(index_mc, mcf = lead(mc) - mc)
index_mc <- mutate(index_mc, mcf = lead(mcf, base - 1))
index_mc <- mutate(index_mc, mcf = lag(mcf, base, 0))
index_mc <- mutate(index_mc, mcf = cumsum(mcf))
index_mc <- arrange(index_mc, -time)
index_mc <- mutate(index_mc, mcb = lead(mc) - mc)
index_mc <- mutate(index_mc, mcb = lead(mcb, nyears - base))
index_mc <- mutate(index_mc, mcb = lag(mcb, nyears - base + 1, 0))
index_mc <- mutate(index_mc, mcb = cumsum(mcb))
index_mc <- group_by(index_mc, time)
index <- summarise(
index_mc,
mean = exp(mean(mcf)) * exp(mean(mcb)),
sd = sd(mcf) * exp(mean(mcf)) + sd(mcb) * exp(mean(mcb))
)
index <- mutate(index, lower = mean - sd, upper = mean + sd)
ggplot(index) +
aes(x = parse_date_time(time, "Y"), y = mean, ymin = lower, ymax = upper) +
geom_ribbon(alpha = .2) +
geom_line() +
ylab(NULL) +
xlab(NULL) +
theme_minimal()