---
output: 
  html_document:
    css: help.css

params:
  compartment: biota
  detGroup: Imposex
---

```{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)

periodStart <- 2010
periodEnd <- 2015
assessmentYear <- 2017

source("report functions.R")

load(file.path("meta analysis objects", "regional data.RData"))
assessment <- get(params$compartment)
assessment <- assessment[[params$detGroup]]
purpose <- assessment$purpose

isAC <- !is.null(assessment$AC)

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

if (isAC) {
  detAClabel <- switch(
    detAC, 
    "Effects Range Low" = "ERL",
    "European Commission food standard" = "EC",
    "Environmental Assessment Criterion" = "EAC"
  )
}

detLabel <- switch(
  params$detGroup, 
  Metals = "metal", 
  Chlorobiphenyls = "CB",
  "PAH (parent)" = "PAH",
  "Organo-bromines" = "BDE",
  "Organo-metals" = "organo-metal",
  Imposex = "VDS"
)


detFixedEffect <- params$detGroup %in% "Metals"

if (params$compartment == "biota" && detFixedEffect) {
  greenID <- "orange" 
  isOrange <- FALSE
} else {
  greenID <- "green"
  isOrange <- any(assessment$status$colour %in% "orange")
}
  
isBlack <- any(assessment$status$colour %in% "black")

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


```

## MIME regional assessment of status and trends of imposex (VDS) in whelks

### Overview

<p class = "startbig">
The `r assessmentYear` MIME assessment of contaminant and biological effects data includes an assessment of the trends and status of imposex (VDS) in whelks at monitoring stations throughout the OSPAR area. Assessments are made for a large number of time series, each of VDS levels in a single species at a single monitoring station. This document synthesises the results of the individual time series to assess the status and trends of VDS at the MIME regional level.  For simplicity, the term 'region' is used throughout to describe MIME regions.  OSPAR regions are always referred to as such.
</p>

<p class = "bot0">
The following sections:
</p>

<ul type = "disc">
<li>summarise the individual time series results</li>
<li>describe the methods used to combine the individual time series results to estimate regional trends and status</li>
<li>present the estimates of regional trends</li>
<li>present the estimates of regional status</li>
</ul>

<br>

### Summary of individual time series results {.tabset}

<p class = "startsmall">
A time series of `r detLabel` concentrations is assessed for status if:
</p>

<ul class = "last">
<li>there is at least one year with data in the period `r periodStart` to `r periodEnd`</li>
<li>there are at least three years of data over the whole time series</li>
</ul>

<p class = "botsmall">
A time series is also assessed for trends if:
</p>

<ul class = "lastbig">
<li>there are at least three years of data with individual measurements or four years of data with pooled measurements</li>
</ul>

<p class = "botbig">
The first tab below shows all the monitoring stations where there are individual time series assessments of trend or status for VDS. The magenta circles are stations where there is a trend assessment, and the cyan cirles are stations where there is only a status assessment.
</p>

<p class = "botsmall">
The next tab show the trend and status assessments for VDS. The colours have the following meaning:
</p>

<ul class = "last">
```{asis, echo = isAC}
<li>blue: mean VDS is significantly (p < 0.05) below the Background Assessment Concentration (BAC)</li>
<li>`r greenID`: mean VDS is significantly (p < 0.05) below the `r detAC` (`r detAClabel`)</li>
<li>red: mean VDS is not signficantly below the `r detAC` (`r detAClabel`)</li>
```

```{asis, echo = isBlack}
<li>black: there are no assessment criteria available</li>
```
</ul>

<p class = "botsmall">
And the shapes have the following meaning:
</p>

<ul type = "disc">
<li>downward triangle: mean VDS is significantly decreasing (p < 0.05)</li>
<li>circle: there is no significant change in mean VDS (p > 0.05)</li>
<li class = "lastbig">upward triangle: mean VDS is significantly increasing (p < 0.05)</li>
</ul>

<p class = "botsmall">
The final set of tabs give:
</p>

<ul class = "lastbig">
<li>the number of time series with upwards, downwards or no trend in each MiME and OSPAR region</li>
<li>the same information presented as proportions</li>
```{asis, echo = isAC}
<li>the number of time series with blue, `r greenID`, or red status in each MIME and OSPAR region</li>
<li>the same information presented as proportions</li>
```
</ul>


<br>

#### Stations

```{r, echo = FALSE, message = FALSE, warning = FALSE, fig.width = 7, fig.height = 7}
plotStations(assessment, purpose = purpose, includeAll = TRUE)
```

<br>


```{r include = FALSE}
out <- NULL
for (detID in assessment$detID) {
  out <- c(out, knit_child('meta_analysis_bydet.Rmd'))
}
```


```{r, echo=FALSE, results="asis"}
cat(paste(out, collapse = '\n'))
```


#### Trend summary (numbers)

<p class = "startbig">
Number of time series with upwards, downwards or no trend by MIME and OSPAR region   
</p>

```{r echo = FALSE, message = FALSE, results = "asis"}
tab <- tableTimeSeries(assessment, "trend", asNumbers = TRUE, purpose = purpose)
tab <- xtable::xtable(tab)
print(tab, type = "html", include.rownames = FALSE) 
```

<br>

#### Trend summary (proportions)

<p class = "startbig">
Proportion of time series with upwards, downwards or no trend by MIME and OSPAR region  
</p>

```{r echo = FALSE, message = FALSE, results = "asis"}
tab <- tableTimeSeries(assessment, "trend", asNumbers = FALSE, purpose = purpose)
tab <- xtable::xtable(tab)
print(tab, type = "html", include.rownames = FALSE) 
```

<br>

```{r, eval = isAC, child = "meta_analysis_text_status_tables.Rmd"}
```


### MIME regional assessment - methods 

<p class = "startbig">
Tabulating the number of time series with each status category by region provides a quick summary of the individual time series results.  However, it does not provide an objective regional assessment of status.  Similarly, tabulating the number of time series with an upward or downward trend does not provide an objective regional assessment of trend.  This section describes how the individual time series results can by synthesised in a meta-analysis to assess both status and trend at the regional level.  
</p>


#### Trends

<p class = "startsmall">
For a regional trend assessment, the trend in each time series is summarised by the estimated odds ratio of the VDS of an individual whelk being above the `r detAClabel` in one year relative to the previous year.  Values of 1 indicate no trend; i.e. the odds of an individual being above the `r detAClabel` in one year are the same as in the year before.  Values < 1 indicate that the odds of being above the `r detAClabel` in one year are lower than in the year before, so there is a decline in the level of imposex.  Conversely, values > 1 indicate an increase in the level of imposex.  Regional trends are then estimated by fitting the following linear mixed model by restricted maximum likelihood:
</p>

<ul class = "last">
<li>response: trend (odds ratio)</li>
<li>fixed model: region</li>
<li>random model: trend estimation variation + residual variation</li>
</ul>

<p class = "botsmall">
The fixed model means that a trend is estimated for each region. The random model has two terms:
</p>

<ul class = "last">
<li>trend estimation variation is the variance of the trend estimates from the individual time series analysis, assumed known and fixed</li>
<li>residual variation is the variation that cannot be explained by any of the fixed effects or the other random effects</li>
</ul>

<p class = "botlarge">
The meta-analysis is restricted to time series from monitoring stations that are classified as 'representative' or 'impacted' (those close to a point source). 'Baseline' stations (those with near pristine conditions or only subject to very diffuse inputs) are omitted because any trends there will likely be caused by different processes to those at representative or impacted stations. Note that the inclusion of impacted stations differs from the approach used for contaminants (where they are excluded). Many imposex monitoring stations were chosen because they were close to harbours (i.e. impacted). However, the ban of TBT (which causes imposex) came into effect at the same time throughout the OSPAR area, so it seems plausible that trends in VDS might be similar (on the logistic scale) in both impacted and representative stations.
</p>

<p class = "botbig">
The meta-analysis is further restricted to regions with at least three trend stations with good geographic spread. Three stations is considered the minimum required to provide some sort of evidence base at the regional level.  
</p>


#### Status

<p class = "startsmall">
For a regional status assessment, the status of each time series is summarised by the ratio between the estimated mean VDS in the final monitoring year and the assessment critrion. This ensures that status is always measured on the same scale, even though the assessment criterion might vary between time series (because VDS is measured in a different species). The ratio is transformed to the square root scale (to better satisfy modelling distributional assumptions) and essentially the same linear mixed model as for trends is then fitted:
</p>


<ul class = "last">
<li>response: status (sqrt(mean VDS / assessment criterion))</li>
<li>fixed model: region</li>
<li>random model: status estimation variation + residual variation</li>
</ul>

<p class = "botbig">
where status estimation variation is the variation in the status estimates from the individual time series analysis, assumed known and fixed.
</p>

<p class = "botbig">
There are no restrictions on the time series used in the status meta-analysis based on the classification of the monitoring station; time series from baseline, representative and impacted stations are all included. Again, the meta-analysis is restricted to regions with at least three status stations with good geographic spread. 
</p>


#### Presentation

<p class = "startbig">
The results of the status regional assessment are back-transformed for presentation. A value of 1 indicates that the regional mean VDS is equal to the assessment criterion; a value of 0.5 indicates that it is half the assessment criterion; and a value of 2 indicates that it is twice the assessment criterion. There is no need to back-transform the estimated regional trend, which is on the odds ratio scale.    
</p>

<br>


### Regional assessment - availability of time series {.tabset}

<p class = "startbig">
This section provides more detail on the number and geographic spread of the time series available for the meta-analyses.  The tabs show the stations with parametric trend and status assessments for each OSPAR region in turn, having excluded trend assessments at baseline and impacted stations.  The magenta circles are stations where there is a trend assessment, and the cyan cirles are stations where there is only a status assessment.
</p>


<p class = "botsmall">

```{r include = FALSE, echo = FALSE}
dropRegionsTrend <- sort(with(assessment, setdiff(levels(trend$region), levels(regionalTrend$region))))
dropRegionsStatus <- sort(with(assessment, setdiff(levels(status$region), levels(regionalStatus$region))))
```

```{asis, echo = length(dropRegionsTrend) > 1 | length(dropRegionsStatus) > 1}
<p class = "botsmall">
The following regions have an insufficient number or time series or spread of stations and are excluded from the regional assessment:
</p>

<ul type = "disc">
<li>trend: `r paste(dropRegionsTrend, collapse = ", ")`</li>
<li>status: `r paste(dropRegionsStatus, collapse = ", ")`</li>
</ul>
```

<br>

```{r include = FALSE}
out <- NULL
for (regionID in as.character(1:4)) {
  out <- c(out, knit_child('meta_analysis_regional_plots.Rmd'))
}
```


```{r, echo=FALSE, results="asis"}
cat(paste(out, collapse = '\n'))
```


### Regional trends {.tabset}

<p class = "startsmall">
The first 2 tabs show:
</p>

<ul class = "last">
<li>Regional trend: regional trend estimates with pointwise 95% confidence limits</li>
<li>Individual time series: the trend estimates from the individual time series on which the regional assessment is based.  The estimates from regions with insufficient geographic spread (and excluded from the meta-analysis) are also shown.</li>
</ul>

<p class = "botsmall">
The symbols in all these plots have the following interpretation:
</p>

<ul class = "lastbig">
<li>downward triangle: the mean VDS is significantly decreasing (p < 0.05)</li>
<li>circle: there is no change in mean VDS (p > 0.05)</li>
<li>upward triangle: the mean VDS is significantly increasing (p < 0.05)</li>
</ul>

```{r include = FALSE, echo = FALSE}
txt <- if(detFixedEffect) "" else
  paste0("regional trend (averaged over ", detLabel, "s) and the ")
```

<p class = "startsmall">
The final tab shows the estimates of the regional trends.  
</p>

<br>

#### Regional trends 

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

<br> 



#### Individual time series

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

<br>

#### Estimates by region 

<p class = "startsmall">
This table shows the estimates of the regional trend:  
</p>

<ul type = "disc">
<li>trend: the estimated mean odds ratio of the VDS of an individual whelk being above the `r detAClabel` in one year relative to the previous year</li>
<li>se: the corresponding standard error</li>
<li>lower, upper: the corresponding pointwise 95% confidence limits</li>
</ul>

<br>

```{r echo = FALSE, message = FALSE, results = "asis"}
tab <- tableTrendImposex(assessment, purpose = purpose, type = "region", detLabel = detLabel)
tab <- xtable::xtable(tab)
print(tab, type = "html", include.rownames = FALSE) 
```

<br>


### Regional status compared to the `r detAC` (`r detAClabel`) {.tabset}

<p class = "startsmall">
The first 2 tabs show:
</p>

<ul class = "last">
<li>Regional status: regional status estimates with pointwise upper one-sided 95% confidence limits</li>
<li>Individual time series: the status estimates from the individual time series on which the regional assessment is based.  The estimates from regions with insufficient geographic spread (and excluded from the meta-analysis) are also shown.</li>
</ul>

<p class = "botsmall">
The symbols in all the plots have the following interpretation:
</p>

<ul class = "last">
<li>blue: the mean VDS is significantly (p < 0.05) below the Background Assessment Concentration (BAC)</li>
<li>green: the mean VDS is significantly (p < 0.05) below the `r detAClabel`</li>
<li>red: the estimate of the mean VDS is not significantly below the `r detAClabel`</li>
</ul>

<p class = "botbig">
Note that to see why some points are coloured blue, it is necessary to look ahead to the next section, where status is compared to the BAC.  Unfortunately, it isn't easy to present the assessment of status relative to the `r detAClabel` with the assessment of status relative to the BAC because the ratio of the BAC to the `r detAClabel` can vary between species.
</p>

<p class = "startsmall">
The final tab shows the estimates of regional status.  
</p>

<br>


#### Regional status

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

<br>


#### Individual time series

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

<br>


#### Estimates by region 

<p class = "startsmall">
This table shows the estimates of the regional status:  
</p>

<ul type = "disc">
<li>status: the estimated mean sqrt VDS relative to the `r detAClabel`)</li>
<li>se: the corresponding standard error</li>
<li>upper: the corresponding pointwise upper one-sided 95% confidence limit</li>
<li>VDS: the estimated mean VDS relative to the `r detAClabel`</li>
<li>VDS upper: the corresponding pointwise upper one-sided 95% confidence limit</li>
</ul>

<br>

```{r echo = FALSE, message = FALSE, results = "asis"}
tab <- tableStatusImposex(assessment, AC = detAClabel, purpose = purpose, type = "region", detLabel = detLabel)
tab <- xtable::xtable(tab)
print(tab, type = "html", include.rownames = FALSE) 
```

<br>


### Regional status compared to the Background Assessment Concentration (BAC) {.tabset}

<p class = "startsmall">
The first 2 tabs show:
</p>

<ul class = "last">
<li>Regional status: regional status estimates with pointwise upper one-sided 95% confidence limits</li>
<li>Individual time series: the status estimates from the individual time series on which the regional assessment is based.  The estimates from regions with insufficient geographic spread (and excluded from the meta-analysis) are also shown.</li>
</ul>

<p class = "botsmall">
The symbols in all the plots have the following interpretation:
</p>

<ul class = "last">
<li>blue: the mean VDS is significantly (p < 0.05) below the Background Assessment Concentration (BAC)</li>
<li>green: the mean VDS is significantly (p < 0.05) below the `r detAClabel`</li>
<li>red: the estimate of the mean VDS is not significantly below the `r detAClabel`</li>
</ul>

<p class = "startsmall">
The final tab shows the estimates of regional status.  
</p>


<br>

#### Regional status

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

<br>



#### Individual time series

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

<br>


#### Estimates by region 

<p class = "startsmall">
This table shows the estimates of the regional status:  
</p>

<ul type = "disc">
<li>status: the estimated mean sqrt VDS relative to the BAC</li>
<li>se: the corresponding standard error</li>
<li>upper: the corresponding pointwise upper one-sided 95% confidence limit</li>
<li>VDS: the estimated mean VDS relative to the BAC</li>
<li>VDS upper: the corresponding pointwise upper one-sided 95% confidence limit</li>
</ul>

<br>


```{r echo = FALSE, message = FALSE, results = "asis"}
tab <- tableStatusImposex(assessment, AC = "BAC", purpose = purpose, type = "region", detLabel = detLabel)
tab <- xtable::xtable(tab)
print(tab, type = "html", include.rownames = FALSE) 
```

<br>


