README.Rmd

---
title: "Modelling Capture Probability"
date: "26 August 2020"
output:
  github_document:
    html_preview: false
---

[![Build Status](https://travis-ci.org/Faskally/ef.svg?branch=tmb)](https://travis-ci.org/Faskally/ef)

## installation

To install the latest version `ef` package run
```r
remotes::install_github("faskally/ef")
```


```{r include=FALSE}
library(ef)
library(mgcv)
library(dplyr)
library(reshape)
library(lattice)
data(ef_data)
```

# Create Required Fields for Analysis


Create a factor for pass with separate levels for the first and subsequent passes. NOTE: Assuming the second and third pass have the same capture probability avoids problems with model identifiability ([Millar _et al_., 2016](https://www.sciencedirect.com/science/article/pii/S0165783616300017))

```{r}
ef_data$pass12 <- factor(replace(ef_data$pass, ef_data$pass > 2, 2))
```

Create a sampleID for grouping EF counts together. A sample is a unique combination of site visit x species x lifestage and is required by the EF function

```{r}
ef_data$sampleID <-
  as.numeric(
    factor(
      paste(ef_data$date, ef_data$siteID, ef_data$lifestage, ef_data$species)
    )
  )
```

Create a visitID for unique site visits. This is required by the overdispersion function (see below).

```{r}
ef_data$visitID <- factor(paste(ef_data$date, ef_data$siteID))
```

Create a factor for year
```{r}
ef_data$fyear <- factor(ef_data$year)
```

# Fit a Model


Fit a simple model. This contains a factor for pass (2 levels), lifestage (2 levels), site (3 levels) and year (5 levels). With a larger dataset it would also be possible to add continuous variables as linear or smoothed terms providing the degrees of freedom are specified (typically up to k=3 (2 d.f.), see [Millar _et al_., 2016](https://www.sciencedirect.com/science/article/pii/S0165783616300017) for further details)

```{r results = 'hide', warning = FALSE, message = FALSE}
m1 <-
  efp(
    count ~ siteID + fyear + pass12 + lifestage,
    data = ef_data, pass = pass, id = sampleID
  )
```
```{r echo = FALSE}
m1
```


# Get Predicted Values for _P_


The model needs to be reformatted to allow predictions (i.e. obtain _P_) capture probabilities, uses as.gam function.

```{r results = 'hide', warning = FALSE, message = FALSE}
m1_gam <- as.gam(m1)
```

Predict _P_ from the model to the full dataset. Note, if all the data is multi-pass then can just use fitted values. However, if you want P for both single and multi-pass data then predict to the new dataframe

```{r}
ef_data$prob <- predict.gam(m1_gam, newdata=ef_data, type = "response")
```


# Obtain Cumulative Estimates of _P_ (for estimation of densities later)


Aggregate counts by sampleID (i.e. counts of fry or parr per site visit)

```{r}
counts <- aggregate(count ~ sampleID, ef_data, sum)
```

Reshape to get capture probability by species x lifestage x site visit x pass. This facilitates the summary of _P_.

```{r results = 'hide', warning = FALSE, message = FALSE}
probs <- cast(ef_data, formula = sampleID ~ pass, value = "prob")
```

Get cumulative capture probability by sampleID (see [Glover _et al_., 2018](https://www.sciencedirect.com/science/article/pii/S1470160X18303534))

```{r}
probs $ prob <- apply(probs[,2:4], 1, function(x) 1 - prod(1-x))
```

Remove individual pass capture probabilities

```{r}
probs <- select(probs, sampleID, prob)
```

Remove pass and pass-specific data from dataframe

```{r}
densdata <- unique(select(ef_data, -count, -prob, -pass, -pass12))
```

Join data and cumulative counts by sampleID

```{r results = 'hide', warning = FALSE, message = FALSE}
densdata <- left_join(densdata, counts, by = "sampleID")
```

Join data and cumulative capture probability by sampleID

```{r results = 'hide', warning = FALSE, message = FALSE}
densdata <- left_join(densdata, probs, by = "sampleID")
```

# Estimate Density



Density (N/m^-2^) is estimated as $\frac{counts}{area*P}$

```{r}
densdata $ Density_Estimate <- densdata $ count / (densdata $ area *
                                                      densdata $ prob)
```

If you want to model counts directly (Poisson model), you need to have cumulative P and area in the offset. For a Poisson model the offset is $\log({area*cumulativeP})$.

```{r}
densdata $ offset <- log(densdata $ area * densdata $ prob)
```

Plot to check

```{r echo = F, fig.width = 6, fig.height = 4.5}
densdata <- densdata[order(densdata $ year), ]

xyplot(Density_Estimate ~ fyear | lifestage, data = densdata, groups = siteID,
       type = 'b', auto.key = T, scales = list(alternating = F),
       xlab = "Year", ylab = expression(paste("Fish m"^-2)))
```


# Comparing Models


Because electrofishing count data are typically overdispersed and because the modelling framework cannot incorpororate random effects, it it necessary to estimate the amount of overdispersion in the data and to account for this when comparing models.

Overdispersion is estimated by fitting a large model that captures most of the systematic variation in the data and comparing this to a visit-wise and saturated model (see [Millar _et al_., 2016](https://www.sciencedirect.com/science/article/pii/S0165783616300017)):

```{r}
large_model <- efp(count ~ fyear + pass12 + lifestage,
                   data = ef_data, pass = pass, id = sampleID)
```

Use the overdispersion function in the ef package and specify the data, visit ID, site ID, sample ID and the large model.

```{r}
od_estimate <- overdispersion(data = ef_data, visitID = "visitID",
                              sampleID = "sampleID",
                              control = list(maxit = 100),
                              largemodel = large_model)
```

View overdispersion estimate. The 'disp' column contains the estimates of within and between sample overdispersion. The largest of these values is used to adjust the BIC.

```{r}
od_estimate
```

Compare models using the adjusted BIC function (Equation 4, [Millar _et al_., 2016](https://www.sciencedirect.com/science/article/pii/S0165783616300017))

  * Full model

```{r}
mfull <- efp(count ~ siteID + year + lifestage + pass12,
             data = ef_data, pass = pass, id = sampleID)

BICadj(mfull, od_estimate)
```

  * Model without lifestage

```{r}
mlife <- efp(count ~ siteID + year + pass12,
             data = ef_data, pass = pass, id = sampleID)

BICadj(mlife, od_estimate)
```

## contributing

After making changes to the code, run:
```r
devtools::document()
```
to update help files and the NAMESPACE

then
```r
devtools::check()
```
to check for errors, and fix warnings and notes

Finally, the readme can be rebuilt using
```r
rmarkdown::render("README.Rmd")
```

then changes can be commited and pushed.