Trend_analysis.Rmd

# Trend Analysis

```{r setup,  echo = F, message=F}

#Load packages
library(knitr)
library(rmarkdown)

```

## Contributor name(s)
Charles Perretti, Geret DePiper, Sean Hardison

## Data name
Time series trend analysis


### Indicator category
Method
<!-- 1. Database pull -->
<!-- 2. Database pull with analysis -->
<!-- 3. Synthesis of published information -->
<!-- 4. Extensive analysis, not yet published -->
<!-- 5. Published methods -->

##Methods
In some cases, time series indicators in the State of the Ecosystem reports are presented with trend lines, although not all time series were selected for trend analysis. As has been shown elsewhere [see @Nicholson2004; @Wagner2013; @VonStorch1999a], power of statistical tests for trend is hampered by low sample size and autocorrelated observations, which are common attributes of SOE indicators. 

In a preliminary simulation study, we explored the consequences of time series length and autocorrelation strength on statistical power of three tests for trend. Our results showed that rates of Type I error, or the false rejection of the null hypothesis when no trend exists, increased for time series at lower sample sizes. A similar result was seen for Type II error rates at low sample sizes (i.e. failure to reject the null hypothesis when trend exists), and both types of error were inflated by autocorrelation. Based on these findings, we selected a minimum series length of N = 30 for indicator time series before assessing for trend. We also chose to use a GLS model selection (GLS-MS) approach to assess for trend in State of the Ecosystem reports, as this approach was least sensitive to Type II error at N &geq; 30. GLS-MS also allowed for both linear and quadratic model fits. 

The model selection process fit four different GLS models to the data using the R package *nlme* [@Pinheiro2017]. The first two models were first-order linear; one was specified with Gaussian error structure, and the second with $AR(1)$ correlation structure. The second two GLS models were quadratic, one specifying Gaussian error and the second an $AR(1)$ correlation structure. All four model fits were compared using second-order AIC (AICc) with the R package *AICcmodavg* [@Mazerolle2017a], and tested against the null hypothesis of no trend through likelihood ratio tests (p < 0.05) against corresponding null models. In SOE time series figures, significant positive trends were colored orange, and negative trends purple. 

### Data source(s)
NA

### Data extraction
NA

### Data analysis


```{r packages, echo=T, message=FALSE, warning=FALSE, include=T}

#R packages
library(dplyr)
library(nlme)
library(AICcmodavg)
library(data.table)
```

<!--Include accompanying R code, pseudocode, flow of scripts, and/or link to location of code used in analyses.-->
```{r analysis, echo = T, message= F, warning=F, include=T}

data.dir <- "./data"
load(file.path(data.dir, "SOE_data_2018.Rdata"))

#--------------------------------GLS Model Selection-----------------------------#
fit_lm <- function(dat) {

  constant_norm <-
    nlme::gls(series ~ 1, 
              data = dat)
  
  constant_ar1 <-
    try(nlme::gls(series ~ 1,
                  data = dat,
                  correlation = nlme::corAR1(form = ~time)))
  if (class(constant_ar1) == "try-error"){
    return(best_lm <- data.frame(model = NA,
                                 aicc  = NA,
                                 coefs..Intercept = NA,
                                 coefs.time = NA,
                                 coefs.time2 = NA,
                                 pval = NA)) 
  } 
  
  
  
  # Linear model with normal error
  linear_norm <- 
    nlme::gls(series ~ time, 
              data = dat)
  
  # Linear model with AR1 error
  linear_ar1 <- 
    try(nlme::gls(series ~ time, 
                  data = dat,
                  correlation = nlme::corAR1(form = ~time)))
  if (class(linear_ar1) == "try-error"){
    return(best_lm <- data.frame(model = NA,
                                 aicc  = NA,
                                 coefs..Intercept = NA,
                                 coefs.time = NA,
                                 coefs.time2 = NA,
                                 pval = NA))
    
  }
  
  # Polynomial model with normal error
  dat$time2 <- dat$time^2
  poly_norm <- 
    nlme::gls(series ~ time + time2, 
              data = dat)
  
  # Polynomial model with AR1 error
  poly_ar1 <- 
    try(nlme::gls(series ~ time + time2, 
                  data = dat,
                  correlation = nlme::corAR1(form = ~time)))
  if (class(poly_ar1) == "try-error"){
    return(best_lm <- data.frame(model = NA,
                                 aicc  = NA,
                                 coefs..Intercept = NA,
                                 coefs.time = NA,
                                 coefs.time2 = NA,
                                 pval = NA))
    
  }
  
  # Calculate AICs for all models
  df_aicc <-
    data.frame(model = c("poly_norm",
                         "poly_ar1",
                         "linear_norm",
                         "linear_ar1"),
               aicc  = c(AICc(poly_norm),
                         AICc(poly_ar1),
                         AICc(linear_norm),
                         AICc(linear_ar1)),
               coefs = rbind(coef(poly_norm),
                             coef(poly_ar1),
                             c(coef(linear_norm), NA),
                             c(coef(linear_ar1),  NA)),
               # Calculate overall signifiance (need to use
               # ML not REML for this)
               pval = c(anova(update(constant_norm, method = "ML"),
                              update(poly_norm, method = "ML"))$`p-value`[2],
                        anova(update(constant_ar1, method = "ML"),
                              update(poly_ar1, method = "ML"))$`p-value`[2],
                        anova(update(constant_norm, method = "ML"),
                              update(linear_norm, method = "ML"))$`p-value`[2],
                        anova(update(constant_ar1, method = "ML"),
                              update(linear_ar1, method = "ML"))$`p-value`[2]))
  
  best_lm <-
    df_aicc %>%
    dplyr::filter(aicc == min(aicc))
  
  
  if (best_lm$model == "poly_norm") {
    model <- poly_norm
  } else if (best_lm$model == "poly_ar1") {
    model <- poly_ar1
  } else if (best_lm$model == "linear_norm") {
    model <- linear_norm
  } else if (best_lm$model == "linear_ar1") {
    model <- linear_ar1
  }
  
  return(list(p = best_lm$pval,
              model = model))
}

#-------------------------------------Plotting code------------------------------------#
soe.plot <- function(data, x.var, y.var, x.label = '', y.label = '', tol = 0.1,
                     x.start = NA, x.end = NA, end.start = 2008, bg.col = background, mean_line = T,
                     end.col = recent, stacked = NA, x.line = 2.5, y.line = 3.5, scale.axis = 1,
                     rel.y.num = 1.5, rel.y.text = 1.5, suppressAxis = FALSE,status  = F,anomaly = F,
                     endshade = TRUE, full.trend = TRUE, point.cex = 1.5, lwd = 2, ymax = TRUE,ymin = TRUE,
                     y.upper = y.upper, y.lower = y.lower, extra = FALSE, x.var2 = x.var2, y.var2 = y.var2,
                     line.forward = FALSE, mean_line.2 = T, cex.stacked = 1, website = T) {
  

  #Select Data
  x <- data[Var == y.var, ]
  x <- x[order(x[, get(x.var)]), ]
  setnames(x, x.var, 'X')
  
  #Set common time step if necessary
  if(is.na(x.start)) x.start <- min(x[, X])
  if(is.na(x.end))   x.end   <- max(x[, X])
  x <- x[X >= x.start, ]
  
  #Set up plot parameters
  if (ymax == TRUE){
    y.max <- max(x[, Value]) + tol * max(x[, Value])
  } else {
    y.max <- as.numeric(y.upper)
  }
  
  if (ymin == TRUE){
    y.min <- min(x[, Value]) - tol * abs(min(x[, Value]))
  } else if (ymin == FALSE){
    y.min <- as.numeric(y.lower)
  }
  
  y.mean <- mean(x[, Value])
  y.sd <- sd(x[, Value])
  
  #Plot blank plot
  plot(x[X >= x.start, list(X, Var)], xlim = c(x.start, x.end),
       ylim = c(y.min,y.max), xlab = '', ylab = '', axes = F, ty = 'n')


  #Add background
  u <- par('usr')
  rect(u[1], u[3], u[2], u[4], border = NA, col = bg.col)
  
  #Add end period shading
  if (endshade == TRUE){
    rect(end.start - 0.5, u[3], u[2], u[4], border = NA, col = end.col)
  }
  
  #Add mean line
  if (anomaly == F){
      if (mean_line == TRUE){
      abline(h = y.mean, col = 'grey', lwd = 3, lty = 2)
      } 
  } else if (anomaly == TRUE){
      abline(h = 0, col = 'grey', lwd = 3, lty = 2)
  }
  
  #Add x y lines
  abline(h = u[3], lwd=3)
  abline(v = u[1], lwd=3)
  
  #Add data points/lines
  points(x[, list(X, Value)], pch = 16, cex = point.cex)
  lines( x[, list(X, Value)], lwd = lwd)
  
  #extra lines
  if (extra == TRUE){
    x2 <- data[Var == y.var2, ]
    x2 <- x2[order(x2[, get(x.var2)]), ]
    setnames(x2, x.var2, 'X2')
    x2 <- x2[X2 >= x.start, ]
    if (mean_line.2 == TRUE){
     abline(h = mean(x2[, Value]), col = 'lightcoral', lwd = 3, lty = 2) 
    }
    points(x2[, list(X2, Value)], pch = 16, cex = point.cex, col = "indianred")
    lines( x2[, list(X2, Value)], lwd = lwd, col = "indianred")
    }
    
  
  #Add axis
  if (suppressAxis == FALSE){
    if(is.na(stacked)) axis(1, cex.axis = 1)
    if(!is.na(stacked)){
      if(stacked!= 'A') axis(3, cex.axis = 1.5, tck = 0.1, labels = F)
    }
  }

  #Stacked axes with 0 overlap so need to remove
  labels <- round((axTicks(2) / scale.axis), 5)
  if(labels[1] == 0) labels[1] <- ''
  axis(2, at = axTicks(2), labels = labels, cex.axis = rel.y.num,
       las = T)

    #Add axis labels
    if(!website){
      if(!is.na(stacked)) text(u[1], u[4], labels = stacked, cex = cex.stacked, adj = c(-0.5, 1.5))
    } else if (website){
      text(u[1], u[4], labels = "", cex = cex.stacked, adj = c(-0.5, 1.5))
    }
    if(is.na(stacked)){
      mtext(1, text = x.label, line = x.line, cex = 1)
      mtext(2, text = y.label, line = y.line, cex = rel.y.text)
    }
  
    if (full.trend == T){
    #Split data into past decade and full time series
    dat <- as.data.frame(x[, list(X, Value)])
    
    dat <- dat %>% dplyr::rename(series = Value) %>%
      mutate(time = seq(1,nrow(dat),1))
    
    # Fit linear model
    lm_out <- fit_lm(dat = dat)
    p <- lm_out$p
    if (p < .05){
        
      newtime <- seq(min(dat$time), max(dat$time), length.out=length(dat$time))
      newdata <- data.frame(time = newtime,
                      time2 = newtime^2)
      lm_pred <- AICcmodavg::predictSE(lm_out$model, 
                                 newdata = newdata,
                                 se.fit = TRUE)

      year <- seq(x$X[1],x$X[length(x$X)],length.out = length(dat$time))

      # Make plot
      if (lm_pred$fit[length(lm_pred$fit)] > lm_pred$fit[1]){
        lines(year, lm_pred$fit, col = main.pos, lwd = 7)
        points(x[, list(X, Value)], pch = 16, cex = point.cex)
        lines( x[, list(X, Value)], lwd = lwd)

        if (line.forward == TRUE){
           lines(year, lm_pred$fit, col = main.pos, lwd = 7)
        }
      } else if (lm_pred$fit[length(lm_pred$fit)] < lm_pred$fit[1]){
        lines(year, lm_pred$fit, col = main.neg, lwd = 7)
        points(x[, list(X, Value)], pch = 16, cex = point.cex)
        lines( x[, list(X, Value)], lwd = lwd)
        if (line.forward == TRUE){
           lines(year, lm_pred$fit, col = main.neg, lwd = 7)
        }
      }
    }
    
    if (extra == TRUE){
      
      # Second variable
      dat <- as.data.frame(x2[, list(X2, Value)])
    
      dat <- dat %>% dplyr::rename(series = Value) %>%
      mutate(time = seq(1,nrow(dat),1))
    
     # Fit linear model
      lm_out <- fit_lm(dat = dat)
      p <- lm_out$p
      points(x2[, list(X2, Value)], pch = 16, cex = point.cex, col = "indianred")
      lines( x2[, list(X2, Value)], lwd = lwd, col = "indianred")
      if (p < .05){
    
        newtime <- seq(min(dat$time), max(dat$time), length.out=length(dat$time))
        newdata <- data.frame(time = newtime,
                      time2 = newtime^2)
        lm_pred <- AICcmodavg::predictSE(lm_out$model, 
                                 newdata = newdata,
                                 se.fit = TRUE)

        year <- seq(x2$X2[1],x2$X2[length(x2$X2)],length.out =length(dat$time))
   
    # Make plot
        if (lm_pred$fit[length(lm_pred$fit)] > lm_pred$fit[1] ){
          lines(year, lm_pred$fit, col = main.pos, lwd = 7)
          points(x2[, list(X2, Value)], pch = 16, cex = point.cex, col = "indianred")
          lines( x2[, list(X2, Value)], lwd = lwd, col = "indianred")
        } else if (lm_pred$fit[length(lm_pred$fit)] < lm_pred$fit[1]){
          lines(year, lm_pred$fit, col = main.neg, lwd = 7)
          points(x2[, list(X2, Value)], pch = 16, cex = point.cex, col = "indianred")
          lines( x2[, list(X2, Value)], lwd = lwd, col = "indianred")
        } 
     }
    }

  }


 
}  


#Add axis labels for stacked plots
soe.stacked.axis <- function(x.label, y.label, x.line = 2.5,rel.x.text = 1.5,
                             y.line = 3.5, rel.y.text = 1.5, outer = TRUE){
  axis(1, cex.axis = rel.x.text)
  mtext(1, text = x.label, line = x.line, cex = rel.x.text, outer = outer)
  mtext(2, text = y.label, line = y.line, cex = rel.y.text, outer = outer)
  
}


#Background colors
background   <- 'white'
recent       <- '#E6E6E6'
main.pos <- rgb(253/255, 184/255, 99/255,  alpha = .9)
main.neg <- rgb(178/255, 171/255, 210/255, alpha = .9)

```

**Example plot**
```{r example figure, echo = T, eval = T}
opar <- par(mar = c(4, 6, 2, 6))

soe.plot(SOE.data.2018, "Time", "Recreational Seafood MA", scale.axis = 10^6,
         end.start = 2008, x.label = 'Year', rel.y.text = 1.5, rel.y.num = 1.1,
         y.line = 2.5, y.label = expression("Fish caught, 10"^6 *" n"))
```



### Further metadata

#### Public availability statement
NA

#### Point of contact
Sean Hardison, sean.hardison@noaa.gov

#### Data steward
NA

### References
<!--List references here-->