---
title: "Using intrinsic and contextual information associated with automated signal detections to improve call recogniser performance"
author: "Nick Leseberg"
date: "`r format(Sys.Date(), '%d %B %Y')`"
output:  
  pdf_document: 
    fig_caption: yes
    fig_height: 4.75
    fig_width: 6.5
    includes:
      in_header: header.tex
    number_sections: yes
    toc: yes
    toc_depth: 4
    keep_tex: true
  word_document: default
  html_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(comment = "",
                      echo = FALSE,
                      message = FALSE,
                      warning = FALSE,
                      out.width = "0.8\\textwidth",
                      fig.align = "center", 
                      fig.height = 5, 
                      fig.width = 7.5)
suppressPackageStartupMessages({
  library(WWRUtilities)
  library(knitr)
  library(tidyverse)
  library(lme4)
  library(PRROC)
  library(ggplot2)
})
options(knitr.kable.NA = '')     ## use a '' for missing values in kable()
to_pdf <- identical(knitr::opts_knit$get("rmarkdown.pandoc.to"), "latex")  ## is this pdf?

getRE <- function(modl) {
  form <- formula(modl) %>%
    deparse() %>%
    paste(collapse = "") %>%
    gsub("[ \t]", "", .) %>% 
    gsub("^[^|]*\\|", "", .) %>%
    gsub("[1() |]", "", .) %>%
    gsub("\\+", " + ", .)
  if(form == "") form <- "fixed effects only"
  return(form)
}

getFE <- function(modl) {
  form <- formula(modl) %>%
    deparse() %>%
    paste(collapse = "") %>%
    sub("^.*~", "", .) %>%
    sub("\\(1 \\|.*$", "", .) %>%
    sub("[+] *$", "", .) %>%
    gsub(" ", "", .) %>%
    gsub("\\+", " + ", .)
  if(form == "1") form <- "intercept only"
  return(form)
}

results <- read.csv("recogniser_results/combined/new_analysis_levels.csv", row.names = 1) %>%
  select(template, call_class, score, file, class, call, period, site, ARU_type, machine, amp_diff) %>%
  data.frame()

```

# Abstract

1. Rapid expansion in the collection of large acoustic datasets to answer ecological questions has generated a parallel requirement for techniques that streamline analysis of these datasets. In many cases, automated signal recognition algorithms, often termed 'call recognisers', are the only feasible option for doing this. To date, most research has focused on what types of recognisers perform best, and how to train these recognisers to optimise performance. 
2. We demonstrate that once recogniser construction is complete and the data processed, further improvements are possible using intrinsic and contextual information associated with each detection. We initially construct a call recogniser for the Night Parrot (Pezoporus occidentalis) using the R package monitoR, and scan a test dataset. We then examine a number of intrinsic variables associated with each detection generated by the recogniser, and several contextual variables, associated with the species' environment and ecology to determine if they might help predict whether a given detection is a true positive (target signal) or false positive (non-target signal). We test several logistic regression models incorporating different combinations of intrinsic and contextual variables, selecting the best-performing model for application. We train the model, using it to calculate the probability each detection is a true or false positive. 
3. Substituting this model-derived probability for raw recogniser score improved the recogniser's performance, reducing the number of detections requiring proofing by 60% to achieve recall of 90%, and by 76% to achieve recall of 75%. 
4. This technique is applicable to any recogniser output, regardless of the underlying algorithm. Application requires an understanding of how the recogniser algorithm determines matches, and knowledge of a species' ecology and environment. Because advanced programming skills and expertise are not required to apply this technique, it will be particularly relevant to field ecologists for whom building and operating call recognisers is an element of their research toolbox, but not necessarily a focus.
 

# Preliminary recogniser assessment 

## Recogniser construction and the assessment dataset

A Night Parrot call recogniser was constructed using the `R` package `MonitoR`, in accordance with procedures outlined in Katz *et al* 2016. To test recogniser performance,  $90\times10 -$minute files containing Night Parrot calls were extracted from an existing database of sound files. The 90 files consisted of 30 files each from known long-term stable roost sites at Mount Edward Graves (`site_1`), Hills (`site_2`), and White Hills (`site_3`). The dataset was also balanced by time, with 10 files from each site taken from `dusk`, `night` and `dawn` periods. `Dusk` files were clipped from the period within one hour of sunset, and commenced at the first call detected within that period. `Dawn` files were clipped from the period within one hour of sunrise, and finished with the final call detected within that period. `Night` files were clipped from any time between the defined `dusk` and `dawn` periods. 

Each 10-minute file was listened to manually, and all Night Parrot calls noted ($n = `r with(results, sum(call))`$), including the `calltype` (per descriptions in Leseberg *et al*. 2019), and whether the call could be detected by visually scanning a spectrogram, or only aurally, by listening to the file (`v_a`). Each 10-minute file was then scanned using the recogniser. Each results file was matched against the manually scanned dataset, with each 'detection' returned by the recogniser (*n* = 31437) assigned a classification (`class`) as either a true positive (`TP`) or false positive (`FP`). Of the 1850 calls detected manually, 110 were not detected by the recogniser. These were added to the dataset and labelled as false negatives (`FN`). Other parameters associated with each detection included: `call` (whether the dataset entry represented a Night Parrot call); `template` (which of the 31 recogniser templates was responsible for the detection);`call_class` (the parent call type for each template); `score` (similarity score assigned by the recogniser); and, `file` (which of the 90 files the detection was extracted from). Added to this dataset were the following parameters for each detection: `period`, `site`, `ARU_type` (whether file was recorded using an `SM3` or `SM4`), `machine` (a unique identifier for the specific ARU used to record the file), and `temp_diff` (for a 1.1 second clip of each detection, the difference between the maximum amplitude and the mean amplitude within the frequency bounds of the template that triggered the detection). This resulted in the following dataframe: 

```{r data_input}
set.seed(20191105)
some_results <- results[sample(nrow(results), 10), ]
some_results <- within(some_results, {
  file <- as.character(file)
  file <- paste( "...", substring(file, nchar(file)-6))
  call <- ifelse(call, "T", "F")
}) %>% arrange(file)

kable(some_results, align = "r", digits = 2, booktabs = TRUE,
      caption = "A selection of data records (slightly abbreviated)")

```

## Exploratory investigation

### Raw performance and the effect of `period` and `site`

In accordance with the recommendations in Knight *et al*. (2017), raw performance of the recogniser was assessed by calculating the area under the curve (AUC) for the resulting precision-recall curve. To examine whether recogniser performance differed according to the factors `period` or `site`, AUC was also calculated separately for those categories. 

```{r roc, fig.height = 5.5, fig.width = 5.5, out.width="0.495\\textwidth", fig.show="hold"}
par(las = 1)

results_pr <- results %>%  
  select(score, call, class, period, site) %>%
  arrange(desc(score)) %>%
  data.frame()

### Calculate data for combined PR curve
results_combine <- results_pr %>%
  within({
    TPs <- cumsum(class == "TP")
    FPs <- cumsum(class == "FP")
    r_score <- TPs/sum(call)
    p_score <- TPs/(TPs+FPs)
    TPs <- NULL
    FPs <- NULL
  })
extra <- within(results_combine[nrow(results_combine), ], {
  p_score <- 0
})
results_combine <- rbind(results_combine, extra)


### Calculate and plot data for period PR curves

results_dusk <- subset(results_pr, period == "dusk") %>%
  within({
    TPs <- cumsum(class == "TP")
    FPs <- cumsum(class == "FP")
    r_score <- TPs/sum(call)
    p_score <- TPs/(TPs+FPs)
    TPs <- NULL
    FPs <- NULL
  })

extra <- within(results_dusk[nrow(results_dusk), ], {
  p_score <- 0
})
results_dusk <- rbind(results_dusk, extra)


results_night <- subset(results_pr, period == "night") %>%
  within({
    TPs <- cumsum(class == "TP")
    FPs <- cumsum(class == "FP")
    r_score <- TPs/sum(call)
    p_score <- TPs/(TPs+FPs)
    TPs <- NULL
    FPs <- NULL
  })

extra <- within(results_night[nrow(results_night), ], {
  p_score <- 0
})
results_night <- rbind(results_night, extra)


results_dawn <- subset(results_pr, period == "dawn") %>%
  within({
    TPs <- cumsum(class == "TP")
    FPs <- cumsum(class == "FP")
    r_score <- TPs/sum(call)
    p_score <- TPs/(TPs+FPs)
    TPs <- NULL
    FPs <- NULL
  })

extra <- within(results_dawn[nrow(results_dawn), ], {
  p_score <- 0
})
results_dawn <- rbind(results_dawn, extra)


make_dashed_square <- local({
  a <- c(0,0,0,1,1)
  b <- c(1,1,0,0,1)
  c <- c(1,0,1,1,0)
  d <- c(0,0,0,1,1)
  function(u = 1, lty = "dashed", col = "grey") {
    segments(u*a, u*b, u*c, u*d, lty = lty, col = col)
    invisible(NULL)
  }
})

setHook("plot.new",
        list(las = function() par(las = 1),
             pch = function() par(pch = 16)),
        "append")

prcurve_periods <- rbind(cbind(time = "combined", results_combine),
                         cbind(time = "dusk",     results_dusk),
                         cbind(time = "night",    results_night),
                         cbind(time = "dawn",     results_dawn)) %>%
  select(period = time, p_score, r_score)

times <- c("combined", "dusk", "night", "dawn")
cols <- c("black", "steel blue", "fire brick", "orange")

aucf <- function(x) with(x, {
  r <- c(0, r_score)
  sum(diff(r)*p_score)
})
areas_raw <- sapply(split(prcurve_periods, prcurve_periods$period), aucf) %>%
  round(4) %>%
  format(nsmall = 4)

plot(p_score ~ r_score, subset(prcurve_periods, period == "combined"),
     xlim=0:1, ylim = 0:1, col = cols[1], xlab = "Recall", ylab = "Precision",
     axes = FALSE, type = "n")
axis(1, at = 0:5/5)
axis(2, at = 0:5/5)
box()
grid()
make_dashed_square()
lines(p_score ~ r_score, subset(prcurve_periods, period == "combined"), col = cols[1])
lines(p_score ~ r_score, subset(prcurve_periods, period == "dusk"), col = cols[2])
lines(p_score ~ r_score, subset(prcurve_periods, period == "night"), col = cols[3])
lines(p_score ~ r_score, subset(prcurve_periods, period == "dawn"), col = cols[4])

oldPar <- par(family = "mono", cex = 0.9)
legend(0.05, 0.22, legend = paste0(format(times), " (", format(areas_raw), ")"),
       lty = c("solid", rep("solid", 2), "solid"), bty="n", text.font = 2,
       col = c("black", "steel blue", "fire brick", "orange"), lwd = 2)
par(oldPar)

### Calculate and plot data for site PR curves

results_site_1 <- subset(results_pr, site == "site_1") %>%
  within({
    TPs <- cumsum(class == "TP")
    FPs <- cumsum(class == "FP")
    r_score <- TPs/sum(call)
    p_score <- TPs/(TPs+FPs)
    TPs <- NULL
    FPs <- NULL
  })

extra <- within(results_site_1[nrow(results_site_1), ], {
  p_score <- 0
})
results_site_1 <- rbind(results_site_1, extra)


results_site_2 <- subset(results_pr, site == "site_2") %>%
  within({
    TPs <- cumsum(class == "TP")
    FPs <- cumsum(class == "FP")
    r_score <- TPs/sum(call)
    p_score <- TPs/(TPs+FPs)
    TPs <- NULL
    FPs <- NULL
  })

extra <- within(results_site_2[nrow(results_site_2), ], {
  p_score <- 0
})
results_site_2 <- rbind(results_site_2, extra)


results_site_3 <- subset(results_pr, site == "site_3") %>%
  within({
    TPs <- cumsum(class == "TP")
    FPs <- cumsum(class == "FP")
    r_score <- TPs/sum(call)
    p_score <- TPs/(TPs+FPs)
    TPs <- NULL
    FPs <- NULL
  })

extra <- within(results_site_3[nrow(results_site_3), ], {
  p_score <- 0
})
results_site_3 <- rbind(results_site_3, extra)

### Plot PR curves

setHook("plot.new",
        list(las = function() par(las = 1),
             pch = function() par(pch = 16)),
        "append")

prcurve_sites <- rbind(cbind(location = "combined",      results_combine),
                       cbind(location = "site_1",         results_site_1),
                       cbind(location = "site_2",       results_site_2),
                       cbind(location = "site_3",  results_site_3)) %>%
  select(location, p_score, r_score)

locations <- c("combined", "site_1", "site_2", "site_3")
cols <- c("black", "steel blue", "fire brick", "orange")

aucf <- function(x) with(x, {
  r <- c(0, r_score)
  sum(diff(r)*p_score)
})
areas_site <- sapply(split(prcurve_sites, prcurve_sites$location), aucf) %>%
  round(4) %>%
  format(nsmall = 4)

plot(p_score ~ r_score, subset(prcurve_sites, location == "combined"),
     xlim=0:1, ylim = 0:1, col = cols[1], xlab = "Recall", ylab = "Precision",
     axes = FALSE, type = "n")
axis(1, at = 0:5/5)
axis(2, at = 0:5/5)
box()
grid()
make_dashed_square()
lines(p_score ~ r_score,
      subset(prcurve_sites, location == "combined"), col = cols[1])
lines(p_score ~ r_score,
      subset(prcurve_sites, location == "site_1"), col = cols[2])
lines(p_score ~ r_score,
      subset(prcurve_sites, location == "site_2"), col = cols[3])
lines(p_score ~ r_score,
      subset(prcurve_sites, location == "site_3"), col = cols[4])

oldPar <- par(family = "mono", cex = 0.9)
legend(0.05, 0.22, legend = paste0(format(locations), " (", format(areas_site), ")"),
       lty = c("solid", rep("solid", 2), "solid"), bty="n", text.font = 2,
       col = c("black", "steel blue", "fire brick", "orange"), lwd = 2)
par(oldPar)

```

The baseline AUC for the recogniser is `r areas_raw[1]`. The precision-recall plots for `period` show that AUC differs markedly between periods. AUC is higher than baseline for clips during the `night` (AUC = `r areas_raw[3]`) suggesting that recogniser performance is improved during this period. AUC is lowest during the `dusk` period (AUC = `r areas_raw[2]`. Improved performance during `night` is likely explained by the lack of other species calling during this period, meaning fewer false positives. 

The precision-recall plots for `site` show no obvious differences across sites.

### Effect of template

Another potential source of variation in recogniser performance is the nature of the templates included in the recogniser. Some Night Parrot calls, such as *dink* and *toot* are relatively simple, short, single-note calls. Short, single notes are a common component of the calls of several species that commonly occur in the study area, increasing the likelihood the templates for these calls will also generate false positives. Conversely, several Night Parrot calls such as *two-note whistle* and *trill* calls are  unique in comparison to other species that occur in the study area, meaning their templates are less likely to generate false positives. To investigate this effect a contingency table was created comparing the quantity of true positive and false positive detections associated with each category of `call_class`.
```{r, results="asis", echo=FALSE}
if(to_pdf) cat("\\newline\n")
```

```{r template}
results_no_FN <- results %>%
  subset(class != "FN") %>%
  within({
    call_class <- factor(call_class)
    template <- factor(template)
  }) %>% 
  select(template, call_class, score, file, class, 
         call, period, site, ARU_type, machine, amp_diff) %>%
  data.frame()
with(results_no_FN, table(factor(call, levels = c(TRUE, FALSE)), call_class)) %>% 
  kable(caption = "Success rate for different template categories", booktabs=TRUE)

```

It is clear that a much higher proportion of false positives are generated by categories of `call_class` associated with simple calls, while the categories associated with more unique calls generate relatively fewer false positives.

# A predictive model to improve recogniser performance

## Modelling concept 

It seems likely that several parameters associated with each detection, in addition to `score`, could be used to assess the likelihood that a detection represents a true positive or false positive detection. A model that incorporates the values of those parameters to predict the probability that a detection is a true positive could be used to improve recogniser performance, with this probability potentially substituted for recogniser score. The model response will be binary (`TP` or `FP`), so will be modelled using a logistic regression. 

## Model definition and selection

The aim of the process will be to use model estimates to predict the probability that any given detection from a future dataset is a Night Parrot. Factors that will be known *a priori* for future datasets will be of the most interest, and will be included as fixed effects. Based on our previous data exploration, a likely fixed effect will be an interaction between `score` and `period`. We can see below, that the relationship between a detection's `score`, and the probability that detection is a true positive, is curved. So, `score` will be fitted as a quadratic term. 

```{r precision-score_relationship, fig.height=4, fig.width=4.5}

results_score <- results_combine %>%
  subset(class != "FN") %>%
  select(score, p_score)

plot((results_score$p_score ~ results_score$score), 
     xlab = "Score", ylab = "Precision", cex = 0.3)
box()
grid()
 
oldPar <- par(family = "mono", cex = 0.8)
par(oldPar)


```

Another likely fixed effect will be `amp_diff`. This is calculated as the difference between the maximum amplitude for a 1.1 second clip of each detection (calculated to include the entire call that triggers the detection), within the frequency bounds of the template that triggered the detection. For example, the `NIPA_how_10` template has lower and upper frequency bounds of 1.97kHz and 2.16kHz respectively. For any detection triggered by this template, `amp_diff` will be the difference between the maximum amplitude within this frequency range, and the mean amplitude within this frequency range. This variable is of interest, because its value is not influenced by whether the detection is a Night Parrot call. So, detections that score highly for both `amp_diff` and `score` are likely to be Night Parrots. A detection that scores highly for `amp_diff` but receives only a moderate or low score, is more likely to be a false positive. A plot of `amp_diff` against `score`, categorised by `class` supports the presence of this relationship. As a result, `amp_diff` will be included as a fixed effect, as an interaction with `score`.

```{r score-amp_diff relationship, fig.height=4, fig.width=5}

with(results, plot(score, amp_diff, col=class, xlab = "Score", ylab = "Amplitude difference", cex = 0.3))
legend('topleft', legend = levels(results$class), col = 1:3, cex = 0.8, pch = 16)
box()
grid()
 
oldPar <- par(family = "mono", cex = 0.9)
par(oldPar)

```

Other potential fixed effects include `call_class` (the parent call type for the template) and `ARU_type` (which model ARU was used).

Other factors that may influence model performance, but will not be known *a priori* for future datasets include `site`, `file` and `machine`. However, including estimates of their variance in the model could improve model performance, so these factors will be modelled as random effects. The factor `template` will be known *a priori*, but contains 31 levels and is nested within `call_class`. Its predictive power will likely be limited, but understanding its impact on model estimates will still be important. It will be included as a potential random effect.

We will use all available data to test our set of candidate models, and determine the best model for call classification. To test model performance we will then partition the dataset, using one third of the test files to train our final model, and the remaining files for validation

### Random effects selection

To determine the most appropriate model, we initially tested a series of models including all possible fixed effects, and combinations of possible random effects. These models were then ranked by AIC.
```{r, results="asis", echo=FALSE}
if(to_pdf) cat("\\newline\n")
```

```{r model-fits, cache=TRUE}

re_mod_1 <- glm(call ~ period*poly(score,2)*amp_diff +
                  ARU_type +
                  call_class,
                family = binomial, data = results_no_FN)

re_mod_2 <- glmer(call ~ period*poly(score,2)*amp_diff +
                    ARU_type +
                    call_class +
                    (1|machine),
                  family = binomial, data = results_no_FN, nAGQ = 0)

re_mod_3 <- glmer(call ~ period*poly(score,2)*amp_diff +
                    ARU_type +
                    call_class +
                    (1|machine) +
                    (1|file),
                  family = binomial, data = results_no_FN, nAGQ = 0)

re_mod_4 <- glmer(call ~ period*poly(score,2)*amp_diff +
                    ARU_type +
                    call_class +
                    (1|machine) +
                    (1|template),
                  family = binomial, data = results_no_FN, nAGQ = 0)

re_mod_5 <- glmer(call ~ period*poly(score,2)*amp_diff +
                    ARU_type +
                    call_class +
                    (1|machine) +
                    (1|site),
                  family = binomial, data = results_no_FN, nAGQ = 0)

re_mod_6 <- glmer(call ~ period*poly(score,2)*amp_diff +
                    ARU_type +
                    call_class +
                    (1|machine) +
                    (1|file) +
                    (1|template),
                  family = binomial, data = results_no_FN, nAGQ = 0)

re_mod_7 <- glmer(call ~ period*poly(score,2)*amp_diff +
                    ARU_type +
                    call_class +
                    (1|machine) +
                    (1|file) +
                    (1|site),
                  family = binomial, data = results_no_FN, nAGQ = 0)

re_mod_8 <- glmer(call ~ period*poly(score,2)*amp_diff +
                    ARU_type +
                    call_class +
                    (1|machine) +
                    (1|template) +
                    (1|site),
                  family = binomial, data = results_no_FN, nAGQ = 0)

re_mod_9 <- glmer(call ~ period*poly(score,2)*amp_diff +
                    ARU_type +
                    call_class +
                    (1|machine) +
                    (1|file) +
                    (1|template) +
                    (1|site),
                  family = binomial, data = results_no_FN, nAGQ = 0)

re_mod_10 <- glmer(call ~ period*poly(score,2)*amp_diff +
                     ARU_type +
                     call_class +
                     (1|file),
                   family = binomial, data = results_no_FN, nAGQ = 0)

re_mod_11 <- glmer(call ~ period*poly(score,2)*amp_diff +
                     ARU_type +
                     call_class +
                     (1|file) +
                     (1|template),
                   family = binomial, data = results_no_FN, nAGQ = 0)

re_mod_12 <- glmer(call ~ period*poly(score,2)*amp_diff +
                     ARU_type +
                     call_class +
                     (1|file) +
                     (1|site),
                   family = binomial, data = results_no_FN, nAGQ = 0)

re_mod_13 <- glmer(call ~ period*poly(score,2)*amp_diff +
                     ARU_type +
                     call_class +
                     (1|file) +
                     (1|template) +
                     (1|site),
                   family = binomial, data = results_no_FN, nAGQ = 0)

re_mod_14 <- glmer(call ~ period*poly(score,2)*amp_diff +
                     ARU_type +
                     call_class +
                     (1|template),
                   family = binomial, data = results_no_FN, nAGQ = 0)

re_mod_15 <- glmer(call ~ period*poly(score,2)*amp_diff +
                     ARU_type +
                     call_class +
                     (1|template) +
                     (1|site),
                   family = binomial, data = results_no_FN, nAGQ = 0)

re_mod_16 <- glmer(call ~ period*poly(score,2)*amp_diff +
                     ARU_type +
                     call_class +
                     (1|site),
                   family = binomial, data = results_no_FN, nAGQ = 0)

```

```{r random_effects_selection}

re_models <- list(re_mod_1, re_mod_2, re_mod_3, re_mod_4, 
               re_mod_5, re_mod_6, re_mod_7, re_mod_8,
               re_mod_9, re_mod_10, re_mod_11, re_mod_12, 
               re_mod_13, re_mod_14, re_mod_15, re_mod_16)

re_AIC <- data.frame("Random effects" = sapply(re_models, getRE), 
                     AIC =              sapply(re_models, AIC),
                     BIC =              sapply(re_models, BIC),
                     Deviance =         sapply(re_models, deviance),
                     "log Lik." =       sapply(re_models, logLik),
                     "Residual DF" =    sapply(re_models, df.residual),
                     check.names = FALSE) %>% arrange(AIC)

kable(re_AIC, row.names = FALSE, digits = c(0, 2, 2, 2, 2, 0),
      caption = "Random effect AIC table", booktabs=TRUE)

```

When ranked by AIC, four models standout as having a much lower AIC than the other 12. These same four models also have a lower BIC than the other twelve models. We will examine the variation in each parameter for these four models.

```{r random_effects_selection_continued, results="asis"}

showMe <- function(vc) {
  nam <- paste0("_", gsub(":", " within ", names(vc)), "_")
  sdv <- sapply(vc, attr, "stddev")
  data.frame(Component= nam, `Std dev` = sdv, check.names = FALSE)
}

re_models <- re_models[order(sapply(re_models, AIC))]

for(i in 1:4) {
  vc <- VarCorr(re_models[[i]])
  print(kable(showMe(vc), row.names=FALSE, digits = 4,
              caption = paste0("__Model ", i, ":__ _",
  as.character(re_AIC$`Random effects`)[i], "_\n"),
              booktabs = TRUE))
}

```

For models containing `machine` or `site`, the contribution of these terms is relatively low, and much lower than for the parameters `file` and `template`. Additionally, the parameter `file` is nested within `machine`, which in turn is nested within `site`. Selecting `file` as a random effect will ensure the variation that would be assigned to the parameters `site` and `machine` will be retained within the random effects structure of the model, without the need to specify these parameters. Therefore, we will drop these terms as random effects, retaining only `file` and `template`. 

The variance associated with `template` is probably explained by differences in the frequency of occurence of different call types. Templates are constructed to cover the range of sound frequency variation within each call type. However, given the small size of the Night Parrot population, and the apparent uniformity in call characteristics of individuals, most calls will be detected by a relatively small number of the templates. Some templates will make few true positive detections, increasing the variance associated with this factor.

The variance associated with `file` is likely to be associated with variation in the strength of calls from each file. Some files contain a number of calls made in proximity to the ARU, which receive higher scores and are more likely to be Night Parrots. In other files, most of the calls are distant from the recorder and receive lower scores. `File` is also likely to contain an amount of variation that cannot be explained by measured contextual variables. Potential contributors include environmental variables such as temperature or humidity, and small amounts of wind noise (even though this was largely controlled for). Other variables could include background noise or the presence of non-target species. 

After dropping `machine` and `site`, the random effects structure for the model will be:

> __`call ~ fixed effects + (1|file) + (1|template)`__

### Fixed effects selection

Next we test the model including all fixed effects and our chosen random effects:

> __`call ~ period*score^2*amp_diff + ARU_type + call_class + (1|file) + (1|template)`__

This model results in the following coefficients:

```{r fixed_effects_selection, cache=TRUE}
fe_mod <- glmer(call ~ period*poly(score,2)*amp_diff +
                  ARU_type +
                  call_class +
                  (1|file) +
                  (1|template),
                family = binomial, data = results_no_FN, nAGQ = 0)

```

```{r fixed_effects_results}
kable(summary(fe_mod)$coeff, digits = 3, booktabs=TRUE,
      caption = "Significance of individual fixed effect coefficients")

kable(dropterm(fe_mod), digits = 3, booktabs=TRUE,
      caption = "Marginal significance of the fixed effects only")

```

These results suggest that the three way interaction between `period`, `score` and `amp_diff` is not influencing model performance, but that each of the two way interactions between these variables should be retained. The influence of `ARU_type` is marginal. `Call_class` has an effect on model performance, but this is not consistent across classes. Calls that are short (i.e *2di*, *1di*, *too* and *ddt*) have less influence on the model than calls which are longer - in each case, the estimate is less than twice the standard error. This is probably related to the increased likelihood of templates for short calls to detect false positives. To investigate this more specifically, we will test the effect of two new variables based on call length. The variable `call_length_1` will categorise all detections based on the template that detects the call as either `short` or `long`, while `call_length_2` will categorise all detections based on the template that detects the call as either `short`, `medium`, or `long`.

We will test a series of models including the original model, plus models retaining only the two way interactions between `period`, `score` and `amp_diff`, a fixed effect representing template category (`call_class`, `call_length_1` or `call_length_2`), and with or without `ARU_type`, meaning seven models will be tested. The random effects structure for each model will include `file` and `template`.

```{r test_fixed_effects, cache=TRUE}

fe_mod_1 <- glmer(call ~ period*poly(score, 2)*amp_diff +
                    call_class +
                    ARU_type +
                    (1|file) +
                    (1|template),
                  family = binomial, data = results_no_FN, nAGQ = 0)

fe_mod_2 <- glmer(call ~ period*poly(score, 2) +
                    poly(score, 2)*amp_diff +
                    period*amp_diff +
                    call_class +
                    ARU_type +
                    (1|file) +
                    (1|template),
                  family = binomial, data = results_no_FN, nAGQ = 0)

fe_mod_3 <- glmer(call ~ period*poly(score, 2) +
                    poly(score, 2)*amp_diff +
                    period*amp_diff +
                    call_class +
                    (1|file) +
                    (1|template),
                  family = binomial, data = results_no_FN, nAGQ = 0)

results_no_FN <- results_no_FN %>% 
  within({
    call_length_1 <- ifelse(call_class %in% c("2di", "1di", "ddt", "too", "3nt"), "short", "long") %>%
      factor(levels = c("short", "long"))
  })
      
fe_mod_4 <- glmer(call ~ period*poly(score, 2)*amp_diff +
                    call_length_1 +
                    ARU_type +
                    (1|file) +
                    (1|template),
                  family = binomial, data = results_no_FN, nAGQ = 0)

fe_mod_5 <- glmer(call ~ period*poly(score, 2) +
                    poly(score, 2)*amp_diff +
                    period*amp_diff +
                    call_length_1 +
                    ARU_type +
                    (1|file) +
                    (1|template),
                   family = binomial, data = results_no_FN, nAGQ = 0)

fe_mod_6 <- glmer(call ~ period*poly(score, 2) +
                    poly(score, 2)*amp_diff +
                    period*amp_diff +
                    call_length_1 +
                    (1|file) +
                    (1|template),
                  family = binomial, data = results_no_FN, nAGQ = 0)

results_no_FN <- results_no_FN %>% 
  within({ 
    call_length_2 <- case_when(call_class %in% c("1di", "2di",
                                                 "too", "ddt")  ~ "brief",
                               call_class %in% c("1tr", "3nt",
                                                 "2wh")         ~ "mid",
                               call_class %in% c("4wh", "2tr",
                                                 "how")         ~ "long") %>%
      factor(levels = c("brief", "mid", "long"))
  })

fe_mod_7 <- glmer(call ~ period*poly(score, 2)*amp_diff +
                    call_length_2 +
                    ARU_type +
                    (1|file) +
                    (1|template),
                  family = binomial, data = results_no_FN, nAGQ = 0)


fe_mod_8 <- glmer(call ~ period*poly(score, 2) +
                    poly(score, 2)*amp_diff +
                    period*amp_diff +
                    call_length_2 +
                    ARU_type +
                    (1|file) +
                    (1|template),
                  family = binomial, data = results_no_FN, nAGQ = 0)

fe_mod_9 <- glmer(call ~ period*poly(score, 2) +
                    poly(score, 2)*amp_diff +
                    period*amp_diff +
                    call_length_2 +
                    (1|file) +
                    (1|template),
                  family = binomial, data = results_no_FN, nAGQ = 0)

```

```{r fixed_effects_analysis}
fe_models <- list(fe_mod_1, fe_mod_2, fe_mod_3, fe_mod_4, fe_mod_5, fe_mod_6, fe_mod_7, fe_mod_8, fe_mod_9)
fe_AIC <- data.frame("Fixed effects" = sapply(fe_models, getFE),
                     AIC =             sapply(fe_models, AIC),
                     BIC =             sapply(fe_models, BIC),
                     Deviance =        sapply(fe_models, deviance),
                     `log Lik.` =      sapply(fe_models, logLik),
                     `Residual DF` =   sapply(fe_models, df.residual),
                     check.names = FALSE)%>% arrange(AIC)
best_modl <- fe_models[[order(sapply(fe_models, AIC))[1]]]

kable(fe_AIC, row.names = FALSE, digits = 2, booktabs=TRUE,
      caption = "Fixed effects AIC table")

```

Ranking models by AIC suggests that `call_class` should be classified according to length using three factors, and `ARU_type` should be retained in the model. This model was marginally better than the model which retained `ARU_type`, but classified template length using two factors. The best ranking model by BIC, which penalises more heavily for complexity, also included `ARU_type`, and classified template length using two factors. When ranked by BIC this model outperformed the model with the best AIC ($\Delta$ BIC < 8), whereas $\Delta$ AIC < 1 between the same two models. Given these parameters, we will retain as the final model:

> __`call ~ period*poly(score,2) + poly(score,2)*amp_diff + period*amp_diff + call_length_1 + ARU_type + (1|file) + (1|template)`__

## Model predictions

We now turn to the test data and see what the chosen model can achieve when applied to new data. We partition the data, using one third of the files to train the model, and setting aside two thirds of the data as the test dataset. Using only the results from the test dataset for comparison, we re-plot our precision-recall curves, substituting the predicted probability for `score`. Because we do not know `file` in advance, we predict using the estimate for this random effect obtained from model training. 

```{r test_detections}

### Randomly select one-third of files as training data, balanced by period
set.seed(20190827)

### select 10 dusk files
fls_dusk <- subset(results, period == "dusk")
fls_dusk <- with(fls_dusk, sort(unique(file)))
train_dusk <- sample(fls_dusk, length(fls_dusk)/3)

### select 10 night files
fls_night <- subset(results, period == "night")
fls_night <- with(fls_night, sort(unique(file)))
train_night <- sample(fls_night, length(fls_night)/3)

### select 10 dawn files
fls_dawn <- subset(results, period == "dawn")
fls_dawn <- with(fls_dawn, sort(unique(file)))
train_dawn <- sample(fls_dawn, length(fls_dawn)/3)

### define training dataset based on sampled files
results_train_dusk <- subset(results,
                             file %in% train_dusk)
results_train_night <- subset(results,
                              file %in% train_night)
results_train_dawn <- subset(results,
                             file %in% train_dawn)

results_train <- rbind(results_train_dusk,
                       results_train_night,
                       results_train_dawn)

### define remaining data as test dataset
results_test_dusk <- subset(results,
                            file %in% setdiff(fls_dusk, train_dusk))
results_test_night <- subset(results,
                             file %in% setdiff(fls_night, train_night))
results_test_dawn <- subset(results,
                            file %in% setdiff(fls_dawn, train_dawn))
results_test <- rbind(results_test_dusk,
                      results_test_night,
                      results_test_dawn)

### remove FNs from training and test datasets

results_train_no_FN <- results_train %>%
  subset(class != "FN") %>%
  within({
    call_class <- factor(call_class)
    template <- factor(template)
  })

results_train_no_FN <- results_train_no_FN %>% 
  within({
    call_length_1 <- ifelse(call_class %in% c("2di", "1di", "ddt", "too", "3nt"), "short", "long") %>%
      factor(levels = c("short", "long"))
  })

results_test_no_FN <- results_test %>%
  subset(class != "FN") %>%
  within({
    call_class <- factor(call_class)
    template <- factor(template)
  }) 

results_test_no_FN <- results_test_no_FN %>% 
  within({
    call_length_1 <- ifelse(call_class %in% c("2di", "1di", "ddt", "too", "3nt"), "short", "long") %>%
      factor(levels = c("short", "long"))
  })

### use training dataset to calibrate model

fe_mod_final <- glmer(call ~ period*poly(score, 2) +
                        poly(score, 2)*amp_diff +
                        period*amp_diff +
                        call_length_1 +
                        ARU_type +
                        (1|template) +
                        (1|file),
                      family = binomial, data = results_train_no_FN, nAGQ = 0)

### use test dataset to get predicted scores
results_test_no_FN$pred <- predict(fe_mod_final,
                                   results_test_no_FN,
                                   re.form = ~(1|template))

FNs <- results_test %>%
  subset(class == "FN")

FNs$pred <- NA
FNs$call_length_1 <- NA

results_test <- rbind(results_test_no_FN, FNs)

```

```{r predicted_roc, fig.height = 5.5, fig.width = 5.5, out.width="0.495\\textwidth", fig.show="hold"}
par(las = 1)

results_pr <- results %>%  
  select(score, call, class, period) %>%
  arrange(desc(score)) %>%
  data.frame()

### Calculate data for score-based PR curves

results_combine <- results_pr %>%
  within({
    TPs <- cumsum(class == "TP")
    FPs <- cumsum(class == "FP")
    r_score <- TPs/sum(call)
    p_score <- TPs/(TPs+FPs)
    TPs <- NULL
    FPs <- NULL
  })
extra <- within(results_combine[nrow(results_combine), ], {
  p_score <- 0
})
results_combine <- rbind(results_combine, extra)


results_dusk <- subset(results_pr, period == "dusk") %>%
  within({
    TPs <- cumsum(class == "TP")
    FPs <- cumsum(class == "FP")
    r_score <- TPs/sum(call)
    p_score <- TPs/(TPs+FPs)
    TPs <- NULL
    FPs <- NULL
  })

extra <- within(results_dusk[nrow(results_dusk), ], {
  p_score <- 0
})
results_dusk <- rbind(results_dusk, extra)


results_night <- subset(results_pr, period == "night") %>%
  within({
    TPs <- cumsum(class == "TP")
    FPs <- cumsum(class == "FP")
    r_score <- TPs/sum(call)
    p_score <- TPs/(TPs+FPs)
    TPs <- NULL
    FPs <- NULL
  })

extra <- within(results_night[nrow(results_night), ], {
  p_score <- 0
})
results_night <- rbind(results_night, extra)


results_dawn <- subset(results_pr, period == "dawn") %>%
  within({
    TPs <- cumsum(class == "TP")
    FPs <- cumsum(class == "FP")
    r_score <- TPs/sum(call)
    p_score <- TPs/(TPs+FPs)
    TPs <- NULL
    FPs <- NULL
  })

extra <- within(results_dawn[nrow(results_dawn), ], {
  p_score <- 0
})
results_dawn <- rbind(results_dawn, extra)

make_dashed_square <- local({
  a <- c(0,0,0,1,1)
  b <- c(1,1,0,0,1)
  c <- c(1,0,1,1,0)
  d <- c(0,0,0,1,1)
  function(u = 1, lty = "dashed", col = "grey") {
    segments(u*a, u*b, u*c, u*d, lty = lty, col = col)
    invisible(NULL)
  }
})

setHook("plot.new",
        list(las = function() par(las = 1),
             pch = function() par(pch = 16)),
        "append")

prcurve_periods <- rbind(cbind(time = "combined", results_combine),
                         cbind(time = "dusk",     results_dusk),
                         cbind(time = "night",    results_night),
                         cbind(time = "dawn",     results_dawn)) %>%
  select(period = time, p_score, r_score)

times <- c("combined", "dusk", "night", "dawn")
cols <- c("black", "steel blue", "fire brick", "orange")

aucf <- function(x) with(x, {
  r <- c(0, r_score)
  sum(diff(r)*p_score)
})
areas <- sapply(split(prcurve_periods, prcurve_periods$period), aucf) %>%
  round(4) %>%
  format(nsmall = 4)

plot(p_score ~ r_score, subset(prcurve_periods, period == "combined"),
     xlim=0:1, ylim = 0:1, col = cols[1], xlab = "Recall", ylab = "Precision",
     axes = FALSE, type = "n")
axis(1, at = 0:5/5)
axis(2, at = 0:5/5)
box()
grid()
make_dashed_square()
lines(p_score ~ r_score, subset(prcurve_periods, period == "combined"), col = cols[1])
lines(p_score ~ r_score, subset(prcurve_periods, period == "dusk"), col = cols[2])
lines(p_score ~ r_score, subset(prcurve_periods, period == "night"), col = cols[3])
lines(p_score ~ r_score, subset(prcurve_periods, period == "dawn"), col = cols[4])

oldPar <- par(family = "mono", cex = 0.9)
legend(0.05, 0.22, legend = paste0(format(times), " (", format(areas), ")"),
       lty = c("solid", rep("solid", 2), "solid"), bty="n", text.font = 2,
       col = c("black", "steel blue", "fire brick", "orange"), lwd = 2)
par(oldPar)

### Calculate and plot data for pred-based PR curves

results_pr_pred <- results_test %>%  
  select(pred, call, class, period) %>%
  arrange(desc(pred)) %>%
  data.frame()

results_combine_pred <- results_pr_pred %>%
  within({
    TPs <- cumsum(class == "TP")
    FPs <- cumsum(class == "FP")
    r_score <- TPs/sum(call)
    p_score <- TPs/(TPs+FPs)
    TPs <- NULL
    FPs <- NULL
  })
extra <- within(results_combine_pred[nrow(results_combine_pred), ], {
  p_score <- 0
})
results_combine_pred <- rbind(results_combine_pred, extra)


results_dusk_pred <- subset(results_pr_pred, period == "dusk") %>%
  within({
    TPs <- cumsum(class == "TP")
    FPs <- cumsum(class == "FP")
    r_score <- TPs/sum(call)
    p_score <- TPs/(TPs+FPs)
    TPs <- NULL
    FPs <- NULL
  })

extra <- within(results_dusk_pred[nrow(results_dusk_pred), ], {
  p_score <- 0
})
results_dusk_pred <- rbind(results_dusk_pred, extra)


results_night_pred <- subset(results_pr_pred, period == "night") %>%
  within({
    TPs <- cumsum(class == "TP")
    FPs <- cumsum(class == "FP")
    r_score <- TPs/sum(call)
    p_score <- TPs/(TPs+FPs)
    TPs <- NULL
    FPs <- NULL
  })

extra <- within(results_night_pred[nrow(results_night_pred), ], {
  p_score <- 0
})
results_night_pred <- rbind(results_night_pred, extra)


results_dawn_pred <- subset(results_pr_pred, period == "dawn") %>%
  within({
    TPs <- cumsum(class == "TP")
    FPs <- cumsum(class == "FP")
    r_score <- TPs/sum(call)
    p_score <- TPs/(TPs+FPs)
    TPs <- NULL
    FPs <- NULL
  })

extra <- within(results_dawn_pred[nrow(results_dawn_pred), ], {
  p_score <- 0
})
results_dawn_pred <- rbind(results_dawn_pred, extra)


setHook("plot.new",
        list(las = function() par(las = 1),
             pch = function() par(pch = 16)),
        "append")

prcurve_pred_periods <- rbind(cbind(time = "combined", results_combine_pred),
                              cbind(time = "dusk",     results_dusk_pred),
                              cbind(time = "night",    results_night_pred),
                              cbind(time = "dawn",     results_dawn_pred)) %>%
  select(period = time, p_score, r_score)

times <- c("combined", "dusk", "night", "dawn")
cols <- c("black", "steel blue", "fire brick", "orange")

aucf <- function(x) with(x, {
  r <- c(0, r_score)
  sum(diff(r)*p_score)
})
areas_model <- sapply(split(prcurve_pred_periods, prcurve_pred_periods$period), aucf) %>%
  round(4) %>%
  format(nsmall = 4)

plot(p_score ~ r_score, subset(prcurve_pred_periods, period == "combined"),
     xlim=0:1, ylim = 0:1, col = cols[1], xlab = "Recall", ylab = "Precision",
     axes = FALSE, type = "n")
axis(1, at = 0:5/5)
axis(2, at = 0:5/5)
box()
grid()
make_dashed_square()
lines(p_score ~ r_score, subset(prcurve_pred_periods, period == "combined"),
      col = cols[1])
lines(p_score ~ r_score, subset(prcurve_pred_periods, period == "dusk"),
      col = cols[2])
lines(p_score ~ r_score, subset(prcurve_pred_periods, period == "night"),
      col = cols[3])
lines(p_score ~ r_score, subset(prcurve_pred_periods, period == "dawn"),
      col = cols[4])

oldPar <- par(family = "mono", cex = 0.9)
legend(0.05, 0.22, legend = paste0(format(times), " (", format(areas_model), ")"),
       lty = c("solid", rep("solid", 2), "solid"), bty="n", text.font = 2,
       col = c("black", "steel blue", "fire brick", "orange"), lwd = 2)
par(oldPar)

```

For the combined data, and for each `period` (`dusk`, `night`, and `dawn`), substituting the model-predicted score for raw score results in an increased AUC, meaning the recogniser performance is improved. As expected, given preliminary analysis, this improvement is modest for the `night` period (AUC = `r areas_model[3]` for model predicted and AUC = `r areas_raw[3]` for raw score), and greatest for the `dawn` period (AUC = `r areas_model[4]` for model predicted and AUC = `r areas_raw[4]` for raw score). 

# P(will be detected if available) for isolated deployments 

Isolated ARUs are commonly deployed in patches of suitable habitat to either monitor known long-term stable roost sites, or in an attempt to detect new populations of the species. In these circumstances the probability that a Night Parrot will be detected if available within the acoustic dataset equals the recogniser's recall when calculated for a dataset where calls may occur at any distance from the ARU, as for the assessment dataset used here.

## Score cutoffs for specified values of recall

Because the distribution of raw scores is not affected by `period` (i.e the raw score of any given call is not affected by `period`), the raw score value that achieves a specified recall is applicable across each `period` for the entire dataset.  However, the distribution of model-predicted scores will be affected by `period`, as this is a fixed effect incorporated into the model. The model parameters suggest that the detections based on certain templates, particularly longer calls, receive higher model-predicted scores than those of shorter calls during the `dusk` and `dawn` periods. This would be comparable to having more birds calling closer to the ARUs during these periods. Therefore recall under the model-predicted scores may vary between periods. Accordingly, the model-predicted score cutoff for each recall threshold was calculated separately for each `period`, using only the test dataset, rather than the entire dataset, to avoid overfitting.

```{r extract_recall_values}

### prepare dataset to extract FPs
efficiency <- results_test %>%
  select(call, class, pred, score, file, period) %>%
  arrange(desc(score)) %>%
  data.frame()

### calculate precision and recall values for score-based data
efficiency_score <- efficiency %>%
  within({
    TPs <- cumsum(class == "TP")
    FPs <- cumsum(class == "FP")
    r_score <- TPs/sum(call)
    p_score <- TPs/(TPs+FPs)
  })

### calculate precision and recall values for dusk pred-based data
efficiency_pred_dusk <- results_test %>%
  select(call, class, pred, score, file, period) %>%
  subset(period == "dusk") %>%
  arrange(desc(pred)) %>%
  within({
    TPs <- cumsum(class == "TP")
    FPs <- cumsum(class == "FP")
    r_pred <- TPs/sum(call)
    p_pred <- TPs/(TPs+FPs)
  })

### calculate precision and recall values for night pred-based data
efficiency_pred_night <- results_test %>%
  select(call, class, pred, score, file, period) %>%
  subset(period == "night") %>%
  arrange(desc(pred)) %>%
  within({
    TPs <- cumsum(class == "TP")
    FPs <- cumsum(class == "FP")
    r_pred <- TPs/sum(call)
    p_pred <- TPs/(TPs+FPs)
  })

### calculate precision and recall values for dawn pred-based data
efficiency_pred_dawn <- results_test %>%
  select(call, class, pred, score, file, period) %>%
  subset(period == "dawn") %>%
  arrange(desc(pred)) %>%
  within({
    TPs <- cumsum(class == "TP")
    FPs <- cumsum(class == "FP")
    r_pred <- TPs/sum(call)
    p_pred <- TPs/(TPs+FPs)
  })

### function to extract raw scores associated with specific recall values
recall_score <- function(score_recall)
  with(subset(efficiency_score, r_score <= score_recall), min(score))
increments <- seq(0.5, 0.9, 0.05)
scores <- Vectorize(recall_score)
score_recall_thresholds <- scores(increments)
score_recall_thresholds_output <- as.data.frame(score_recall_thresholds)

### remove TPs for FP calculations
efficiency_score_FPs <- subset(efficiency_score, class == "FP")

### function to extract number of FPs with score > recall threshold
score_detections <- function(score_recall_threshold)
  with(efficiency_score_FPs, tapply(score, file, function(x)
    sum(x > score_recall_threshold)))
score_total_detections <- Vectorize(score_detections)
score_total_detections <- score_total_detections(score_recall_thresholds)

### compile results into dataframe
colnames(score_total_detections) <- paste0(substring(100+increments, 3))
score_total_detections <- cbind(with(results_no_FN, {
  data.frame(file = sort(unique(as.character(file))),
             period = tapply(as.character(period), file, unique)) }),
  score_total_detections)
rownames(score_total_detections) <- NULL
score_total_detections <- score_total_detections %>%
  drop_na()


### function to extract pred scores associated with recall thresholds
### for dusk period
recall_pred <- function(pred_recall)
  with(subset(efficiency_pred_dusk, r_pred <= pred_recall), min(pred))
preds <- Vectorize(recall_pred)
pred_recall_thresholds_dusk <- preds(increments)
pred_recall_thresholds_dusk_output <- as.data.frame(pred_recall_thresholds_dusk)

### remove TPs for FP calculations
efficiency_pred_dusk_FPs <- subset(efficiency_pred_dusk, class == "FP")

### function to extract number of FPs with pred > recall threshold for dusk
pred_detections <- function(pred_recall_threshold)
  with(efficiency_pred_dusk_FPs, tapply(pred, file, function(x)
    sum(x > pred_recall_threshold)))
pred_total_detections_dusk <- Vectorize(pred_detections)
pred_total_detections_dusk <- pred_total_detections_dusk(pred_recall_thresholds_dusk) %>%
as.data.frame() %>%
  drop_na()
colnames(pred_total_detections_dusk) <- paste0(substring(100+increments, 3))

### function to extract pred scores associated with specific recall
### for night period
recall_pred <- function(pred_recall)
  with(subset(efficiency_pred_night, r_pred <= pred_recall), min(pred))
preds <- Vectorize(recall_pred)
pred_recall_thresholds_night <- preds(increments)
pred_recall_thresholds_night_output <- as.data.frame(pred_recall_thresholds_night)

### remove TPs for FP calculations
efficiency_pred_night_FPs <- subset(efficiency_pred_night, class == "FP")

### function to extract number of FPs with pred > recall threshold for night
pred_detections <- function(pred_recall_threshold)
  with(efficiency_pred_night_FPs, tapply(pred, file, function(x)
    sum(x > pred_recall_threshold)))
pred_total_detections_night <- Vectorize(pred_detections)
pred_total_detections_night <- pred_total_detections_night(pred_recall_thresholds_night) %>%
as.data.frame() %>%
  drop_na()
colnames(pred_total_detections_night) <- paste0(substring(100+increments, 3))

### function to extract pred scores associated with specific recall
### for dawn period
recall_pred <- function(pred_recall)
  with(subset(efficiency_pred_dawn, r_pred <= pred_recall), min(pred))
preds <- Vectorize(recall_pred)
pred_recall_thresholds_dawn <- preds(increments)
pred_recall_thresholds_dawn_output <- as.data.frame(pred_recall_thresholds_dawn)

### remove TPs for FP calculations
efficiency_pred_dawn_FPs <- subset(efficiency_pred_dawn, class == "FP")

### function to extract number of FPs with pred > recall threshold for dawn
pred_detections <- function(pred_recall_threshold)
  with(efficiency_pred_dawn_FPs, tapply(pred, file, function(x)
    sum(x > pred_recall_threshold)))
pred_total_detections_dawn <- Vectorize(pred_detections)
pred_total_detections_dawn <- pred_total_detections_dawn(pred_recall_thresholds_dawn) %>%
as.data.frame() %>%
  drop_na()
colnames(pred_total_detections_dawn) <- paste0(substring(100+increments, 3))

### print raw score and model-predicted score cutoffs for recall thresholds
score_cutoffs <- cbind(score_recall_thresholds_output,
                       pred_recall_thresholds_dusk_output,
                       pred_recall_thresholds_night_output,
                       pred_recall_thresholds_dawn_output)
score_cutoffs$recall <- increments
score_cutoffs <- score_cutoffs[,c(5,1,2,3,4)]

cutoffs_table <- kable(score_cutoffs, col.names = c(" ",
                                                    " ",
                                                    "Dusk",
                                                    "Night",
                                                    "Dawn"),
                          format = "latex",
                          row.names = FALSE, digits = 2, booktabs=TRUE,
                          caption = "Raw score and model-predicted score cutoffs for specified recall thresholds")

out_cutoffs <- kableExtra::add_header_above(cutoffs_table, 
                                    c("Recall" = 1,
                                      "Raw Score" = 1,
                                      "Model-derived Probability" = 3))
sub("rrrrr", "@{}ccrrr@{}", out_cutoffs)

```

## Improvements from model application

Recall is related to precision, and this relationship manifests itself in the number of detections that must be checked to achieve the desired level of recall. Using a modelling approach to predict the probability that a detection is a Night Parrot and substituting this for the raw recogniser score significantly improved the performance of the recogniser, improving precision for given values of recall, and reducing the number of detections that must be checked to achieve the desired level of recall.

To understand the practical impact of this improvement, we calculated the mean number of false positive detections that needed to be checked per 10-minute file in the test dataset, to achieve a specific recall. The number of false positive detections is a direct proxy for the amount of time an analyst needs to spend proofing recogniser output. This was calculated by applying the raw score and model-predicted score cutoffs calculated previously for specific values of recall, and calculating the number of false positive detections returned using that cutoff.

Because recogniser performance is known to vary between `dusk`, `night` and `dawn` periods, the mean number of false positive detections per 10-minute clip was determined for each `period`. We also combined these values to determine how many detections would need to be checked to analyse a typical 12-hour night, with one hour of `dusk`, 10 hours of `night`, and 1 hour of `dawn`.

```{r efficiency_analysis}

### get mean number of detections for recall = x and period = dusk
score_total_detections_dusk <- subset(score_total_detections, period == "dusk")
mean_detections_score_dusk <- as.data.frame(sapply(score_total_detections_dusk,
                                                   mean)) %>%
  drop_na()

mean_detections_pred_dusk <- as.data.frame(sapply(pred_total_detections_dusk,
                                                  mean)) %>%
  drop_na()
 
### get mean number of detections for recall = x and period = night
score_total_detections_night <- subset(score_total_detections, period == "night")
mean_detections_score_night <- as.data.frame(sapply(score_total_detections_night,
                                                   mean)) %>%
  drop_na()

mean_detections_pred_night <- as.data.frame(sapply(pred_total_detections_night,
                                                  mean)) %>%
  drop_na()

### get mean number of detections for recall = x and period = dawn
score_total_detections_dawn <- subset(score_total_detections, period == "dawn")
mean_detections_score_dawn <- as.data.frame(sapply(score_total_detections_dawn,
                                                    mean)) %>%
  drop_na()

mean_detections_pred_dawn <- as.data.frame(sapply(pred_total_detections_dawn,
                                                   mean)) %>%
  drop_na()

### combine mean_detections dataframes
mean_detections <- cbind(mean_detections_score_dusk,
                         mean_detections_pred_dusk,
                         mean_detections_score_night, 
                         mean_detections_pred_night, 
                         mean_detections_score_dawn,
                         mean_detections_pred_dawn)

mean_detections$Recall <- increments
mean_detections <- mean_detections[,c(7,1,2,3,4,5,6)]

### calculate number of 'detections' to be checked for a single 12-hour night
### and add to mean_detections dataframe. The total number of 10-minute
### clips will be determined by calculating the total possible calling period
### within the night (with calling potentially starting at sunset, and finishing
### at sunrise. Dusk model will be used within 60 mins after
### sunset, the dawn model within 60 minutes of sunrise, and the
### night model for the remainder
night_length_clips <- ((12*60)-120)/10 
dusk_clips <- 6
dawn_clips <- 60/10
night_clips <- night_length_clips - dusk_clips - dawn_clips

mean_detections$'score' <- mean_detections[,2]*dusk_clips +
  mean_detections[,4]*night_clips + mean_detections[,6]*dawn_clips
mean_detections$'pred' <- mean_detections[,3]*dusk_clips +
  mean_detections[,5]*night_clips + mean_detections[,7]*dawn_clips
mean_detections$'%' <- signif(((mean_detections[,9]/mean_detections[,8])*100),2)

efficiency_table <- kable(mean_detections, col.names = c("Recall",
                                                         "Score", "MDP",
                                                         "Score", "MDP",
                                                         "Score", "MDP",
                                                         "Score", "MDP", "%"),
                          format = "latex",
                          row.names = FALSE, digits = 2, booktabs=TRUE,
                          caption = "Mean false positive detections per 10-minute clip across all periods, and for a simulated 12-hour night")

out <- kableExtra::add_header_above(efficiency_table, 
                                    c(" ",
                                      "Dusk" = 2,
                                      "Night" = 2,
                                      "Dawn" = 2,
                                      "12-hour Night" = 3))
sub("rrrrrrrrrr", "@{}cccccccccc@{}", out)

full_night_score_90 <- round(mean_detections[9,8],0)
full_night_pred_90 <- round(mean_detections[9,9],0)
percent_effort_90 <- mean_detections[9,10]

full_night_score_80 <- round(mean_detections[7,8],0)
full_night_pred_80 <- round(mean_detections[7,9],0)
percent_effort_80 <- mean_detections[7,10]

```

These results clearly demonstrate the efficiencies gained through using model-predicted scores to determine cutoffs for analysis. For a full night of recordings, to achieve recall of 90% using raw score cutoffs would require checking a minimum of `r full_night_score_90` detections. Achieving the same recall using model-predicted cutoffs would require checking `r full_night_pred_90` detections, around `r percent_effort_90`% of the effort if raw score cutoffs were used. Efficiency is even more pronounced if recall of 80% is satisfactory; using model-predicted cutoffs would require approximately `r percent_effort_80`% of the effort if raw score cutoffs were used (`r full_night_score_80` detections for raw score vs `r full_night_pred_80` for model-predicted score). 


