---
title:  "rtrim 2.0 extensions"
author: Patrick Bogaart
date:   "`r Sys.Date()`"
output: 
  rmarkdown::html_vignette:
    toc: true
vignette: >
  %\VignetteIndexEntry{rtrim 2.0 extensions}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
---

```{r, echo = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width  = 7,
  fig.height = 5
)
rm(list=ls())
```

# Introduction

This vignette describes the extensions to the `rtrim` package as per version 2.0:

* Monthly data
* Extended plotting of time-totals
* Extended plotting of indices
* Running TRIM in a stratified manner
* Heatmap visualisations of both observed and imputed data

and in addition extensions introduced in version 2.2:
* long output
* scaled index standard errors

See also the vignettes [rtrim confidence intervals](rtrim_confidence_intervals.html) and [taming overdispersion](taming_overdispersion.html) for additional new features.

# Monthly data

One of the major improvements in `rtrim` 2.0 with respect to 1.0 is the handling of monthly, or any other intra-annual, data.
For example, where a classic TRIM model 3 looks like
$$ \ln \mu_{ij} = \alpha_i + \beta_j $$ 
where $\mu_{ij}$ is an expected count, $\alpha_i$ is a site parameter for site $i$ and $\beta_j$ is a time point parameter for year $j$.
the extension towards months looks like
$$ \ln \mu_{ijm} = \alpha_i + \beta_j + \delta_m $$
where $\delta_m$ are month parameters for month $m$.

Please take a look in [Models and statistical methods in rtrim](TRIM_methods_v2.pdf) for a full explanation, and application to other model formulations.

The general syntax to specify R-Trim models that use monthly data is as follows:
```r
  z <- trim(count ~ site + year, data=...)                   # simple annual data
  z <- trim(count ~ site + year + habitat, data=...)         # using covariates
  z <- trim(count ~ site + (year+month), data=...)           # using monthy data
  z <- trim(count ~ site + (year+month) + habitat, data=...) # using both
```
Note the use of brackets to distinguish months from covariates.

Here is a full example for Oystercatcher data, which now comes with `rtrim` 2.0
```{r}
rm(list=ls()) # always start with a clean slate
library(rtrim)
data(oystercatcher)
oc <- trim(count ~ site + (year+month), data=oystercatcher, model=3, overdisp=TRUE,
           constrain_overdisp=0.999)
plot(index(oc))
```

## Comparison with UIndex

While in the past TRIM was used to analyse count data with an annual resolution (i.e. one observation per site per year), the software package UIndex [Underhill, 1989; Underhill and Prŷs-Jones, 1995] was and is used to analyse count data with higher (e.g., monthly) resolution.
As demonstrated above, `rtrim` is extended to accept and analyse monthly data as well. This section demonstrates the application of `rtrim` to monthly data, and compares the output with that of UIndex.

### UIndex

UIndex was used to analyse the monthly Oystercatcher counts, collected by SOVON Netherlands. Here we show the pre-saved output of UIndex, as the main trend, and the 90\% `consistency intervals'. Note also the use of 2004 as base year.
```{r}
load("UIndex_Oystercatcher_output.RData")
yrange <- range(uidx$index, uidx$lower, uidx$upper)
plot(uidx$year, uidx$index, type='l', xlab="Year", ylab="Index", ylim=yrange)
segments(uidx$year, uidx$lower, uidx$year,uidx$upper)
legend("topright", "UIndex", col="black", lty="solid")

# Add index=1 line for reference
abline(h=1.0, lty="dashed", col=gray(0.5))

# Mark the base year
ibase <- match(2004, uidx$year)
points(uidx$year[ibase], uidx$index[ibase], pch=16)
```

### `rtrim`

Here we show the comparision with `rtrim`, using the results computed above.

```{r}
# Compute and plot an index for Oystercatcher counts, using 2004 as base year and
# adding 90% confidence intervals as well.
idx <- index(oc, level=0.9, base=2004)
plot(idx, band="ci")

# Plot UIndex on top
lines(uidx$year, uidx$index)
segments(uidx$year, uidx$lower, uidx$year,uidx$upper, lwd=2)

legend("bottom", c("UIndex","TRIM"), col=c("black","red"), lty="solid")
```

Note the computation and display of confidence intervals, which is new for `rtrim` 2.0, along with the standard errors of both classic TRIM and `rtrim` 1.0.

This plot demonstrates that the indices as computed by UIndex and `rtrim` are virtually identical, and that the 90% confidence intervals of TRIM are well comparable to the 90% consistency intervals of TRIM, although both are estimated using completely different approaches. In the case of TRIM, confidence intervals are based on standard errors which are derived analytically as part of the GEE estimation process and ultimately are based on the variance within the orginal data.
See the vignettes [Models and statistical methods in rtrim](TRIM_methods_v2.pdf) and [rtrim confidence intervals](rtrim_confidence_intervals.html) for more information. Consistency intervals in UIndex are estimated by means of a bootstrap method. See Underhill [1989] and Underhill and Prŷs-Jones [1995] for more information.

# Stratified `rtrim`

Sometimes it can be usefull to combine `rtrim` results for different regions ('strata') into a single, larger scale ('superstratum') `rtrim` analysis. One particular example is the case where individual EU countries use TRIM or `rtrim` to compute indices for their own countries, and submit the results to the European Bird Census Counsil [EBCC](https://www.ebcc.info/) for aggregatation on the EU scale, see van Strien et al. [2001] for an example using the previous stand-alone version of TRIM.
In this case, the output of the lower scale `rtrim` runs, i.e., the time totals, are used as 'observations' in the higher scale run.
Although this procedure works out well for the estimates and indices, it doesn't work for the associated standard errors, because the time totals are not Poisson distributed, where the original counts are. To circumvent this problem, `rtrim` has options to export the variances of the lower scale runs and to import these into the higher scale runs, to use instead.

The following example shows the associated workflow. Strictly for demonstration purposes, we split the Skylark dataset into two 'regions' associated with the habitats (heath and dunes).

```{r}
# split data
data(skylark2)
heath <- subset(skylark2, habitat=="heath") # 208 records
dunes <- subset(skylark2, habitat=="dunes") # 232 records

heath$site <- factor(heath$site) # get rid of empty levels
dunes$site <- factor(dunes$site)

# run models
m1 <- trim(count ~ site + year, data=heath, model=3)
m2 <- trim(count ~ site + year, data=dunes, model=3)

# collect imputed time-totals (which is the default)
t1 <- totals(m1)
t2 <- totals(m2)

plot(t1,t2, names=c("heath", "dunes"))
```
Note the use of multiple time-totals in a single plot (new for `rtrim` 2.0)

The next step is to use the time totals for the differente habitats (`strata') as  input data for an upscaled (`superstratum') run. The habitat types now serve as site names, and imputed counts will be the input counts.
```{r}
t1$region <- "heath"
t2$region <- "dunes"
t12 <- rbind(t1, t2)
head(t12)
```

The final preparation step is to extract the variance-covariance information for the different habitats, and combine them into a single list, using habitat/region names as identifier, enabling the correct match between the site identifiers in the data, and the variance-covariance matrices.
```{r}
# Also collect the variance-covariance matrices for both runs
vcv1 <- vcov(m1)
vcv2 <- vcov(m2)
vcv3 <- list(heath=vcv1, dunes=vcv2)
```

and off we go with the superstratum run. Note the new argument `covin` to use the variance-covariance data.
```{r}
m3 <- trim(imputed ~ region + time, data=t12, model=3, covin=vcv3)
plot(totals(m3))
```

Now, just for comparison, we compare index plots for both the baseline run (where dunes and heath are taken together, but do act as covariates) and the upscaled `superstratum' variant.
```{r}
m0 <- trim(count ~ site + year + habitat, data=skylark2, model=3) # baseline
t0 <- totals(m0)
t3 <- totals(m3)
plot(t0,t3, names=c("baseline","superstrata"))
```

Which suggests that for *this* example the differences are small, if any.


# Taming overdispersion

In some cases, especially with clustering bird species, overdispersion can be huge, reaching unrealistic values of more than 500.
rtrim now contains an option to constrain the computed value of overdispersion by detecting outliers, and removing them from the computation of overdispersion (but retaining them for all other calculations). The full rationale and methdology is described in [Taming overdispersion](taming_overdispersion.html), but the actual use is rather simple.

Take for example the Oystercatcher data, which results in a huge overdispersion of about 850
```{r}
data(oystercatcher)
m1 <- trim(count ~ site + (year + month), data=oystercatcher, model=3, overdisp=TRUE)
overdispersion(m1)
```

The inclusion of the option `constrain_overdisp=0.999` triggers the detection of outliers that have a probability of 0.1%.
```{r}
m2 <- trim(count ~ site + (year + month), data=oystercatcher, model=3, overdisp=TRUE,
           constrain_overdisp=0.99)
overdispersion(m2)
```

And so we get a much more reasonable result, with smaller standard errors.
```{r}
t1 <- totals(m1)
t2 <- totals(m2)
plot(t1, t2, names=c("unconstrained","constrained"), leg.pos="bottom")
```

# Output visualization

## Plotting time-totals

Once an `rtrim` model has been estimated, one of the first steps of analysis schould be the plotting of time-totals.
This is done by first calling the `totals()` function, and then a custom `plot()` function:

```{r}
rm(list=ls())                          # New section; time for a new blank slate
data(skylark2)                                             # reload Skylark data
m1 <- trim(count ~ site + year, data=skylark2, model=3)
t1 <- totals(m1)          # By default, the time-totals for the imputed data set
plot(t1)
```

Alternatively, one may compute the fitted time-totals. the next example shows the plotting of both the imputed and fitted time-totals, and also demonstrates how series can be named, and the plot can be decorated with a main title.

```{r}
m2 <- trim(count ~ site + year, data=skylark2, model=2, changepoints=c(1,2))
ti <- totals(m2, "imputed")
tf <- totals(m2, "fitted")
plot(ti, tf, names=c("imputed","fitted"), main="Skylark", leg.pos="bottomright")
```

Since imputed totals are composed of both observed and estimated counts, it might be insightful to plot the observed counts as well.
```{r}
m3 <- trim(count ~ site + year, data=skylark2, model=3)
t3 <- totals(m3, obs=TRUE)          # Extract observations in addition to totals
plot(t3)
```

As can be seen, the amount of observed Skylarks is considerable smaller than the time totals suggest.
Furthermore, it can be seen that while the observed counts decrease from 1989, the imputed counts continue to increase.
It is thus suggested to look into more detail what is going on in different sites.


## Plotting indices

Once a TRIM model has been estimated, and indices are computed, these latter can be plotted using the generic plot command `plot()` (which, behind the screens, calls `plot.trim.index()`).

```{r}
m <- trim(count ~ site + year, data=skylark2, model=3) # Run a fairly basic TRIM model
idx <- index(m) # By default, the indices for the imputed data set
plot(idx)
```

If required, the x-axis and y-axis labels as well as the tile can be defined,
and the index can be expressed as a percentage, instead as a fraction.
This example shows all these options:
```{r}
plot(idx, xlab="Year AD", ylab="Index (%)", main="Skylark index", pct=TRUE)
```

### Indices and covariates

When covariates are involved, it can be helpful to compute and plot indices for the various covariate categories as well.
The next example demonstrates this.

```{r}
m <- trim(count ~ site + year + habitat, data=skylark2, model=3) # Run a fairly basic TRIM model
idx <- index(m, covars=TRUE)
plot(idx)
```

As can be seen, indices for the various covariate categories are automatically plotted as well. This behaviour can be supressed by setting `covar="none"` in the call to `plot()`
(note the use of plural `covars' in the call to `index()` --- because indices for multiple covariates can be computed, and the singular `covar` in the call to `plot()` --- because only one of them can be used for a single figure)

### Combining multiple indices

Indices for multiple TRIM runs can be combined in a single plot.

```{r}
data(skylark2)
m0 = trim(count ~ site + year          , data=skylark2, model=3)
m1 = trim(count ~ site + year + habitat, data=skylark2, model=3)

idx0 <- index(m0)
idx1 <- index(m1)

plot(idx0, idx1)
```
As you see, a legend is inserted automatically. You can change the names of the series by using the `names` argument:
```{r}
plot(idx0, idx1, names=c("Without covariates", "Using 'Habitat' as covariate"))
```

## Adding confidence intervals.

New in rtrim 2.0 is the possibility to express uncertainty as a confidence interval, in addition to the standard errors.
Both the functions `totals()` and `index()` now accept the option `level` that specifies the confidence level and triggers the computation.
```{r}
m <- trim(count ~ site + year, data=skylark2, model=3)
tt <- totals(m, level=0.95)                   # Compute 95% confidence intervals
head(tt)
```

So, the lower and upper bounds of the confidence interval is stored in columns `lo` and `hi`. These are automatically picked up by the `plot()` function.
```{r}
plot(tt)
```

If required, the uncertainty band, which is by default plotted using standard errors, can be plotted using the confidence intervals when the option `band="ci"` is used.
```{r}
plot(tt, band="ci")
```

See vignette [rtrim confidence intervals](rtrim_confidence_intervals.html) for more information on the underlying methodology.

## Long output

Stating with version 2.2, the `totals()` and `index()` functions can also generate so-called *long* output, which simplifies plotting using alternative approaches, e.g. `ggplot2` or similar, especially because imputed and fitted time-totals are in rows, not in columns.
```{r}
tt <- totals(m, long=TRUE)
head(tt)

lo <- tt$value - 1.96 * tt$SE # Assume normal distribution
hi <- tt$value + 1.96 * tt$SE

# create an empty plot with sufficient space
xrange = range(tt$year)
yrange = range(lo, hi)
plot(xrange, yrange, type='n', xlab="Year", ylab="Time-totals")
# plot the time series and error bars
segments(tt$year, lo, tt$year, hi, col="red", lwd=2)
lines(tt$year, tt$value, col="red", lwd=2)
```

A similar mechanism has been included for custom plotting of overall trends, using a new separate `trendlines` function.
```{r}
tl <- trendlines(overall(m))     # collect overall trend line
print(tl)

tt <- totals(m, long=TRUE)       # collect time-totals

# define plot limits
xr <- range(tt$year)
yr <- range(tl$lo, tl$hi, tt$value)
plot(xr, yr, type='n', xlab="Year", ylab="Total counts")

# Plot uncertainty band
ubx <- c(tl$year, rev(tl$year))
uby <- c(tl$lo, rev(tl$hi))
polygon(ubx, uby, col=gray(0.9), border=NA)

# Plot trend line
lines(tl$year, tl$value, col="black", lwd=2)

# Plot time-totals
lines(tt$year, tt$value, col="red", lwd=2)
points(tt$year, tt$value, col="red", pch=16, cex=1.5)
```
## "Scaled" standard errors

In `rtrim` version $<2.2$ standard errors for indices are always computed using a formal approach that results in $SE=0$ for the reference year.
Mathematically, this makes completely sense because the indices are computed by dividing the time-totals time series by the time-total for the reference year. The index for that reference year is thus by definition 1, without any uncertainty, hence $SE=0$. The following plot illustrates this.
```{r}
tt  <- totals(m)
idx <- index(m)
par(mfrow=c(1,2))
plot(tt, main="Time-totals", ylab=NA)
plot(idx, main="Index", ylab=NA)
```

However, for many (communication) purposes this formal approach is rather confusing,
because the interpretation is often made that uncertainties 'disappeared' this way,
Furthermore, the indices are often only calculated to compare trends for multiple species.
For this comparison purposes, one often wants to preserve the original uncertainty pattern as in the time-totals.
Starting with `rtrim` version 2.2, this is possible by using the `method="scaled"` option:
```{r}
fidx <- index(m, method="formal") # same as just index(m)
sidx <- index(m, method="scaled")
par(mfrow=c(1,2))
plot(fidx, main="Formal approach")
plot(sidx, main="Scaled approach")
```
As can be seen, the uncertainty for the reference year now is $>0$, while the uncertainties for the other years have been shrunk, which can be explaiend from what loosely can be described as a preservation of total uncertainty.


## Plotting heatmaps.

The detailed spatiotemporal structure of both the observed and the imputed data can be inspected by means of the function `heatmap()` that operates on the output of `trim()`. The default behaviour of this function is to display a heat map of the observed counts only:

```{r}
m <- trim(count ~ site + year, data=skylark2, model=3)
heatmap(m, main="Skylark, observations")
```
In this heatmap, site/time combinations are colored by (log) counts: lower counts are colored a pale red, and high counts a dark red. Consistent with the nature of count data, this color scale is proportional to the log counts. Observed counts of 0 cannot be represented this way and are colored white. Missing site/time combinations are marked as gray.

It can be seen that the observational coverage is not constant: most sites have incomplete records, especially in the earlier years. This is a typical patern for an expanding observation program, but may have consequences for the statistical analysis, because the imputation for these years will in fact be an extrapolation back in time. 

The next example shows the TRIM estimated counts (using shades of blue, rather than red:
```{r}
heatmap(m, "fitted", main="Skylark, TRIM estimates")
```
From this plot, it is clear that the variance between sites is much higher than the variance between years.
In fact, the trend in time can hardly be seen.

The last example sows the heatmap for the imputed data, where estimates are used to fill up the missing observations.
Again, red is for obervations, blue for estimates.
```{r}
heatmap(m, "imputed", main="Skylark, imputed data")
```

### Heatmaps for monthly datasets

For monthly data, heatmaps work slightly different, but in the same spirit:

```{r}
data(oystercatcher)
m <- trim(count ~ site + (year + month), data=oystercatcher, model=3, overdisp=TRUE)
heatmap(m, "imputed", main="Oystercatcher (imputed)")
```

Again, observational coverage is extremely variable in both space and time.
There appears to be a few sites that have sporadic, yet high, count observations, causing large amounts of estimated counts for this location for all other time points, which may effect the aggregated time-totals in a significant way.

Also note that in this example, many site/time combinations have registered a count of 0, which are colored white, as explained above.
<!-- which is possible from an observational point of view, but difficult to reconcile with the log-linear trend assumptions made by TRIM. In this heatmap, these cases are colored white. -->

# References

van Strien, A. J., J. Pannekoek and D.W. Gibbons (2001) Indexing European bird population trends using results of national monitoring schemes: a trial of a new method, Bird Study, 48 (2), 200-213, DOI: [10.1080/00063650109461219](http://doi.org/10.1080/00063650109461219)

Underhill, LG, Prŷs-Jones, RP (1994) Index numbers for waterbird populations. I. Review and methodology. J Appl Ecol, 31, 463-480. doi: [10.2307/2404443](https://www.jstor.org/stable/2404443)

Underhill, L.G. (1989) Indices for Waterbird Populations. BTO Research Report 52, British Trust for Ornithology, Tring.

