---
title: "Appendix 6. Code to develop and analyze noise models for a single species"
author: '.. et al.' 
date: 'March 20, 2020'
output:
  pdf_document:
    fig_height: 3
    fig_width: 5
  html_document:
    fig_height: 3
    fig_width: 5
  word_document:
    fig_height: 3
    fig_width: 5
editor_options: 
  chunk_output_type: console
---


Load packages
```{r, include=FALSE}
# This changes the default colors in lattice plots.
require(mosaic)
require(ggformula)
require(tidyverse)
require(gamm4)
require(MuMIn)
# if you don’t have devtools, instructions to install it are at:
# https://www.r-project.org/nosvn/pandoc/devtools.html
devtools::install_github('stacyderuiter/s245')
require(s245)
knitr::opts_chunk$set(
 tidy=FALSE,     # display code as typed
 size="small"# slightly smaller font for code
)
```

Build myQAIC function
```{r}
myQAIC <- function(model, ...){
  # thanks to Ben Bolker:
  # https://cran.r-project.org/web/packages/bbmle/vignettes/quasi.pdf
  if (!missing(...)){#if multiple models
    dots.xq <- lapply(list(model, ...), function(x) stats::update(x, family='x.quasipoisson'))
    dots.noq <- lapply(list(model, ...), function(x) stats::update(x, family='poisson'))
    chats <- lapply(dots.noq, s245::dfun)
    result <- mapply(MuMIn::QAIC,
                     object=dots.xq,
                     chat = chats)
    Call <- match.call()
    Call$k <- NULL
    names(result) <- as.character(Call[-1L])
    result <- data.frame(result)
      }else{#if only one model
    xq.model <- stats::update(model, family='x.quasipoisson')
    noq.model <- stats::update(model, family='poisson')
    result <- MuMIn::QAIC(xq.model, chat = dfun(noq.model))
  }
  return(result)
}
```

Load dataset
```{r}
allspecies <- read_csv("allspecies_11Dec18.csv")
```

Filter dataset for species of interest
```{r}
#set your species here using the four letter species code (https://www.birdpop.org/docs/misc/Alpha_codes_eng.pdf)
#example here is for the American goldfinch (AMGO)
myspec <- allspecies %>%
  filter(SPEC=='AMGO') %>% # <---------- Change 'YOURSPECIES' to the correct 4-letter code!!
  arrange(YEAR, STATION) %>% dplyr::filter(YEAR != "All")
```


##Productivity models
Prepare the dataset:
```{r}
# function to prep prod dataset
prep_prod_data <- function(data){
  data <- data %>%
    dplyr::select(SPEC, DECLAT, DECLNG, YEAR, productivity, wtmean, sd, STATION, 
                  PC1, PC2, PC3, PC4, PC5, PC6, PC7, PC8, ADULTS) %>%
    mutate(JUVENILES = ADULTS*productivity) %>%
    mutate(SPEC = factor(SPEC),
           STATION = factor(STATION),
           YEAR = as.numeric(as.character(YEAR))) %>%
    na.omit() %>% # this effectively removes the year 'All'
    arrange(YEAR, STATION)
  return(data)
}

myspec_prod <- prep_prod_data(myspec)  
```


Fit the full model:
```{r}
#Count data models (quasi-Poisson) for the number of juveniles with an offset for the number of adults.

prodmod <- gam(JUVENILES ~ wtmean + sd + s(DECLAT, DECLNG, k=20)+ s(YEAR, k=8) + 
                        PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC7 + PC8 + 
                  s(STATION, bs='re') + offset(log(ADULTS)), method='ML', select=TRUE,
                family=quasipoisson(link='log'), data=myspec_prod, na.action='na.fail')

summary(prodmod)
plot(resid(prodmod))
acf(resid(prodmod))
```

Develop and rank candidate models
```{r}
myQAIC(prodmod, update(prodmod, .~.-wtmean), update(prodmod, .~.-sd),
       update(prodmod, .~.-wtmean - sd))
```

Reporting: 
Strength of evidence for noise effect options are: 

- strong, top model includes noise effect, and no other model within 2 AIC units lacks noise effect.
- weak, top model has noise effect but others within 2 AIC units do not
- none, top model has no noise


## Abundance model

Prep dataset
```{r}
#function to prep abun data
prep_abun_data <- function(data){
  data <- data %>%
    dplyr::select(SPEC, DECLAT, DECLNG, YEAR, ADULTS, net.hrs, wtmean, sd, STATION, adults600hr,
                  PC1, PC2, PC3, PC4, PC5, PC6, PC7, PC8)%>%
    dplyr::filter(net.hrs != 0) %>%
    mutate(SPEC = factor(SPEC),
           STATION = factor(STATION),
           YEAR = as.numeric(as.character(YEAR))) %>%
    na.omit() %>%
    arrange(YEAR, STATION)
  return(data)
}

myspec_abun <- prep_abun_data(myspec) 
```

Fit the full model:
```{r}
#Poisson distribution with log link for abundance 

abunmod <- gam(ADULTS ~ wtmean + sd + s(DECLAT, DECLNG, k=20)+ s(YEAR, k=8) + 
                        PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC7 + PC8 + 
                  s(STATION, bs='re') + offset(log(net.hrs)), method='ML', select=TRUE,
                family=quasipoisson(link='log'), data=myspec_abun, na.action='na.fail')

summary(abunmod)
plot(resid(abunmod))
acf(resid(abunmod)) 
```

Develop and rank candidate models
```{r}

myQAIC(abunmod, update(abunmod, .~.-wtmean), update(abunmod, .~.-sd),
       update(abunmod, .~.-wtmean - sd))
```

Reporting:
Strength of evidence for noise effect options are: 

- strong, top model includes noise effect, and no other model within 2 AIC units lacks noise effect.
- weak, top model has noise effect but others within 2 AIC units do not
- none, top model has no noise