---
title: "Single-taxon indicator (rbms)"
output: 
  rmarkdown::html_vignette:
    toc: true
vignette: >
  %\VignetteIndexEntry{Single-taxon indicator (rbms)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---



This documents outlines how a single-taxon indicator (_Aglais urticae_
on farms in Finland) can be calculated in R using the RBMS method.

## Load packages

The following packages are required. All packages are available on
[CRAN](https://cran.r-project.org/) apart from `{fbi}` and `{rbms}` which can be
installed from GitHub.


```r
library(dplyr)
library(fbi)
library(finbif)
library(ggplot2)
library(lubridate)
library(rbms)
```

## Survey data

These five fields are required for the survey data.


```r
select <- c("document_id", "location_id", "year", "month", "day")
```

These filters restrict the survey data to the "Butterflies in Finnish
agricultural landscapes" dataset with the selected data fields have no missing
data.


```r
filter <- list(
  collection = "Butterflies in Finnish agricultural landscapes",
  has_value = select
)
```

The survey data can now downloaded from FinBIF.


```r
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.


```r
surveys <- format_date(surveys)
```

## Count data

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`).


```r
select <- c("document_id", "section", abundance = "abundance_interpreted")
```

The count data requires the same filters as the survey data (though the filter
`has_value` needs to be redefined).


```r
filter[["has_value"]] <- select
```

The count data for _Aglais urticae_ (Small tortoiseshell) can now be downloaded
from FinBIF.


```r
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.


```r
counts <- sum_over_sections(counts, surveys)

counts <- combine_with_surveys(counts, surveys)
```

# Fit RBMS Model

An RBMS model is used to estimate the change in abundance over time (see 
[RBMS documentation
](https://retoschmucki.github.io/rbms/articles/Get_Started_2.html) for details)
with uncertainty estimated via bootstrapping.


```r
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
)
```

# Create Index

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.


```r
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)
```


```r
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()
```

![Relative abundance of Small tortoiseshell on Farms in Finland (calculated using RBMS)](figure/plot-single-taxon-rbms-1.png)
