---
title: "D8: meta-analysis of trends and status of organo-bromine concentrations in biota"
author: "Rob Fryer"
date: "August 1, 2016"
output: html_document
---

```{r setup, include = FALSE, echo = FALSE}

library(knitr)
knit_engines$set(asis = function(options) {
  if (options$echo && options$eval) knit_child(text = options$code)
})


options(width = 110)

compartment <- "biota"
detGroup <- "Organo-bromines"

periodStart <- 2009
periodEnd <- 2014
assessmentYear <- 2015

source("functions.R")

load(file.path("meta analysis objects", "stations by detGroups.RData"))
assessment <- get(compartment)
assessment <- assessment[[detGroup]]

detAC <- switch(
  compartment, 
  sediment = switch(
    detGroup, 
    Metals = "Effects Range Low", 
    Chlorobiphenyls = "Environmental Assessment Criterion",
    "PAH (parent)" = "Effects Range Low",
    "missing"
  ),
  biota = switch(
    detGroup, 
    Metals = "European Commission food standard", 
    Chlorobiphenyls = "Environmental Assessment Criterion",
    "PAH (parent)" = "Environmental Assessment Criterion",
    "missing"
  )
)

detAClabel <- switch(
  detAC, 
  "Effects Range Low" = "ERL",
  "European Commission food standard" = "EC",
  "Environmental Assessment Criterion" = "EAC", 
  "missing" = "CONC"
)

detLabel <- switch(
  detGroup, 
  Metals = "metal", 
  Chlorobiphenyls = "CB",
  "PAH (parent)" = "PAH",
  "Organo-bromines" = "BDE",
  organoMetal = "organo-metal"
)

detUnit <- switch(
  compartment, 
  sediment = switch(
    detGroup, 
    Metals = "mg/kg dw", 
    Chlorobiphenyls = "ug/kg dw",
    "PAH (parent)" = "ug/kg dw",
    "Organo-bromines" = "ug/kg dw",
    organoMetal = "ug/kg dw"
  ),
  biota = switch(
    detGroup,
    Metals = "ug/kg ww", 
    Chlorobiphenyls = "ug/kg ww",
    "PAH (parent)" = "ug/kg ww",
    "Organo-bromines" = "ug/kg ww",
    organoMetal = "ug/kg ww"
  )
)

detFixedEffect <- detGroup %in% "Metals"

isBiota <- compartment == "biota"

textNormalisation1 <- function(compartment, detGroup)
  switch(
    compartment, 
    sediment = paste(switch(detGroup, Metals = "aluminium", "organic carbon"), 
                     "concentrations are also available for normalisation."), 
    biota = "auxiliary data are available to transform concentrations to a wet weight basis if necessary."
  )

textNormalisation2 <- function(compartment, detGroup)
  switch(
    compartment, 
    sediment = paste("normalised to", switch(detGroup, Metals = "5% aluminium", "2.5% organic carbon")),
    biota = "expressed on a wet weight basis"
  )
```

<br>

#### Overview

The `r assessmentYear` assessment of CSEMP data (soon to be available on the BODC website) describes the trends and status of contaminant concentrations in biota and sediment at monitoring stations around the UK.  Assessments are made for a large number of time series, each of a single contaminant in a single species (for biota) at a single monitoring station.  This document is one of a series that synthesises the results of the CSEMP assessment to assess trends and status at the CP2 regional level.  In particular, this document presents a meta-analysis of trends and status of `r detLabel` concentrations in `r compartment`.  The `r detLabel`s considered are `r paste(assessment$detID, collapse = ", ")`.

<br>

#### Trend

A time series of `r detLabel` concentrations is assessed for trends if:

* there are at least five years of data each containing at least one measurement above the detection limit
* there is at least one year with data in the period `r periodStart` to `r periodEnd`

The map below shows all the monitoring stations for which there is a trend assessment for at least one `r detLabel`.  Stations are coloured blue if baseline, green if representative, and red if impacted.  The areas are CP2 regions.

```{r, echo = FALSE, message = FALSE, warning = FALSE, fig.width = 7, fig.height = 7}
plotStations.CSEMP(assessment, "trend")
```

The trend in each time series is summarised by the estimated change in log concentration over the last twenty years (or shorter if the time series doesn't extend that far back). For interpretation, this is back-transformed to the percentage yearly change in concentration.  The plot below shows the estimated precentage yearly changes by `r detLabel` and CP2 region.  All station types are included (baseline, representative and impacted) and coloured accordingly.  

```{r, echo = FALSE, message = FALSE, warning = FALSE, fig.width = 9, fig.height = 5}
plotTrend(assessment)
```

<br>

The number of stations with trend data is tabulated below by region and determinand.  `r if (compartment == "biota") "Note that the number of time series might be more than this, because sometimes several species are monitored at the same station."`

```{r, echo = FALSE, comment = ""}
tableStationsTS(assessment$trend)
```

<br>

The meta-analysis of trends only considers:

* representative stations: impacted stations (those close to a point source) are omitted because they are 'unrepresentative' of contaminant levels at the regional scale; baseline stations (near pristine conditions or only subject to very diffuse inputs) are also omitted because any trends there will likely be caused by different processes to those at representative stations
* CP2 regions with at least three stations with good geographic spread: three stations is considered the minimum required to provide some sort of evidence base at the regional level  
* `r detLabel`s with trend information for at least three stations (over all regions): `r detLabel`s are treated as random effects in the meta-anlaysis, so information about trends in one `r detLabel` can be 'borrowed' from the trend information in the other `r detLabel`s

```{asis, eval = !isBiota}
Unfortunately, this doesn't leave enough data to do anything meaningful with.
```



```{asis, eval = isBiota}
This leaves the following stations:
```

```{r, echo = FALSE, comment = "", eval = isBiota}
tableStationsTS(assessment$regionalTrend)
```


```{asis, eval = isBiota}
<br>

Regional trends are assessed by fitting the following linear mixed model by restricted maximum likelihood:

* response: trend (yearly change in log concentration)
* fixed model: region
* random model: `r detLabel` + region.`r detLabel` + station + trend estimation variation + residual variation

The fixed model means that a trend is estimated for each region.  Trends in the different `r detLabel`s are assumed to be broadly similar, because inputs and control measures are likely to affect all `r detLabel`s in a broadly similar way.  Each regional trend can thus be interpreted as the common trend in `r detLabel` concentrations in that region (i.e. the trend averaged across all `r detLabel`s).  

The random model has five terms:

* `r detLabel` allows for variation in trend between `r detLabel`s common to all regions
* region.`r detLabel` allows for additional variation in trend between `r detLabel`s within regions
* station allows for variation in trend between stations common to all `r detLabel`s
* trend estimation variation is the variation in the trend estimates from the analysis of the original time series, assumed known and fixed
* residual variation is the variation that cannot be explained by any of the fixed effects or the other random effects

<br>

The estimates of regional trends, with pointwise 95% confidence limits, are shown below. The colours signify:

* red: significant upward trend (p < 0.05)
* orange: non-significant trend (p > 0.05)
* green: significant downward trend (p < 0.05)
```


```{r, echo = FALSE, message = FALSE, warning = FALSE, fig.width = 9, fig.height = 5, eval = compartment == "biota"}
plotTrendResults(assessment, "region")
```

```{asis, eval = isBiota}

<br>
Differences in the regional trends between `r detLabel`s is explored below by plotting the estimated trend for each region by `r detLabel`.
```

```{r, echo = FALSE, message = FALSE, warning = FALSE, fig.width = 9, fig.height = 5, eval = isBiota}
plotTrendResults(assessment, "region by determinand")
```

```{asis, eval = isBiota}
Here are the same results but plotted by region
```

```{r, echo = FALSE, message = FALSE, warning = FALSE, fig.width = 9, fig.height = 5, eval = isBiota}
plotTrendResults(assessment, "determinand by region")
```

<br>
<br>

#### Status relative to Assessment Criteria

There are no Assessment Criteria for `r detLabel`s in `r compartment`

<br>
<br>

#### `r detLabel` profiles

The  profile of concentrations across `r detLabel`s can sometimes be useful for examining sources, and for comparing concentrations across regions where Assessment Criteria are not available. The profiles are based on a meta-analysis of the estimated mean log concentration in the final monitoring year from each time series.  The mean log concentration in the final monitoring year is estimated for a time series of `r detLabel` concentrations if:

* there are three or more years of data, at least two of which contain at least one measurement above the detection limit
* there is at least one year with data in the period `r periodStart` to `r periodEnd`

The map below shows all the monitoring stations for which the mean log concentration is estimated for at least one `r detLabel`.  Stations are coloured blue if baseline, green if representative, and red if impacted.  The areas are CP2 regions.

```{r, echo = FALSE, message = FALSE, warning = FALSE, fig.width = 7, fig.height = 7}
plotStations.CSEMP(assessment, "status")
```

The plot below shows the estimated mean log concentrations (backtransformed to the concentration scale) by `r detLabel` and CP2 region.  All station types are included (baseline, representative and impacted) and coloured accordingly.  

```{r, echo = FALSE, message = FALSE, warning = FALSE, fig.width = 9, fig.height = 5}
plotStatus(assessment, detAClabel)
```

<br>

The number of stations with estimated mean log concentration data is tabulated below by region and determinand:

```{r, echo = FALSE, comment = ""}
tableStationsTS(assessment$status, detAClabel)
```

<br>

The concentration profiles are based on a meta-analysis that only considers:

* baseline and representative stations: impacted stations (those close to a point source) are omitted because they are 'unrepresentative' of contaminant levels at the regional scale
* CP2 regions with at least three stations with good geographic spread: three stations is considered the minimum required to provide some sort of evidence base at the regional level  
* `r detLabel`s with estimated mean log concentrations for at least three stations (over all regions): `r detLabel`s are treated as random effects in the meta-anlaysis, so information about the mean log concentration in one `r detLabel` can be 'borrowed' from the estimated mean log concentrations in the other `r detLabel`s

This leaves the following stations:

```{r, echo = FALSE, comment = ""}
tableStationsTS(assessment$regionalStatus, detAClabel)
```

<br>

The concentration profiles are estimated by fitting the following linear mixed model by restricted maximum likelihood:

* response: mean log concentration 
* fixed model: region 
* random model: `r detLabel` + region.`r detLabel` + station + estimation variation + residual variation

The region fixed effect estimates a  mean log concentration for each region 'averaged' across `r detLabel`s, and can be used to compare overall levels across regions.  The `r detLabel` profiles are assumed to be broadly similar across regions, because inputs and control measures are likely to affect all `r detLabel`s in a broadly similar way. This is accommodated by treating `r detLabel` and the interaction between `r detLabel` and region as random effects.  There are five terms in the random model:

* `r detLabel` allows for variation in mean log concentration between `r detLabel`s common to all regions
* region.`r detLabel` allows for additional variation in mean log concentration between `r detLabel`s within regions
* station allows for variation in mean log concentration between stations common to all `r detLabel`s
* estimation variation is the variation in the estimates of mean log concentration from the analysis of the original time series, assumed known and fixed
* residual variation is the variation that cannot be explained by any of the fixed effects or the other random effects

The regional mean concentrations averaged across `r detLabel`s are shown below.  The units are `r detUnit`.

```{r, echo = FALSE, message = FALSE, warning = FALSE, fig.width = 9, fig.height = 5}
plotStatusResults(assessment, "region", "CONC")
```

The estimated `r detLabel` concentration profiles by region are shown below with pointwise 95% confidence limits.  The units are `r detUnit`.


```{r, echo = FALSE, message = FALSE, warning = FALSE, fig.width = 9, fig.height = 5}
plotStatusResults(assessment, "determinand by region", "CONC")
```