---
title: "Variograms and Model Selection"
author: "Christen H. Fleming and Justin M. Calabrese"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Variograms and Model Selection}
  %\VignetteKeyword{variogram}
  %\VignetteEngine{knitr::rmarkdown}
  \usepackage[utf8]{inputenc}
---

In this vignette, we walk through data preparation, variogram analysis, and maximum likelihood estimation.

# Data Preparation

We highly recommend that you get your data onto [Movebank](http://www.movebank.org/).
This will help ensure that your data are of the correct format for `ctmm`, help you identify outliers, and you can keep your data completely private if you wish.
`ctmm` requires that your dataframe conforms to [Movebank naming conventions](http://rstudio.com) (see `help(as.telemetry)`).
The next step is then to import your MoveBank csv file:
```
yourAnimals <- as.telemetry("yourAnimalsMoveBank.csv")
```
Alternatively, if you want to clean your csv file first, you can import it as a data frame
```
yourAnimalsDF <- read.csv("yourAnimalsMoveBank.csv")
```
and then edit the data frame before converting it into a `telemetry` object for `ctmm` via
```
yourAnimals <- as.telemetry(yourAnimalsDF)
```
`as.telemetry` also works on `Move` objects, which can be useful as the `move` package interfaces directly with MoveBank through `R` (see `help(move::getMovebankData)`).

For most species, the default two-point equidistant projection will be fine, however you can provide any PROJ.4 formatted projection with the `projection` argument (see `help(as.telemetry)`). A single fixed projection should be used if you are going to plot groups of individuals that span multiple MoveBank files. This is done by default if multiple individuals are included in a single data frame.

The output of `as.telemetry` will be an individual `telemetry` object or list of `telemetry` objects, depending on how many individual animals are in your csv file. The basic structure of a `telemetry` object is a data frame with columns `t` for time in seconds, and `x` and `y` for the projected locations in meters.

Our example buffalo data is already prepared into a list of `telemetry` objects.
Let us look at the first buffalo and then every buffalo:

```{r,  fig.show='hold'}
library(ctmm)
data("buffalo")
Cilla <- buffalo$Cilla
plot(Cilla)
title("1 Buffalo")
plot(buffalo,col=rainbow(length(buffalo)))
title("5 Buffalo")
```

Looking at the raw movement tracks is a good way to pick out any obvious migratory behaviors.
In the future, we will have migration models to select, but for now all of our models are range resident and so only those portions of the data should be selected.
These buffalo all look fairly range resident, and so we can move on to variograms.

# Variograms

Variograms are an unbiased way to visualize autocorrelation structure when migration, range shifting, drift, or other translations of the mean location are not happening. When drift occurs in the data, then the variogram represents a mixture of both the drift and the autocorrelation structure, each of which contains distinct movement behaviors. In the future, we will have models that can allow for drift, but the current models assume range residence, which we can check with the variogram.

```{r,  fig.show='hold'}
SVF <- variogram(Cilla)
level <- c(0.5,0.95) # 50% and 95% CIs
xlim <- c(0,12 %#% "hour") # 0-12 hour window
plot(SVF,xlim=xlim,level=level)
title("zoomed in")
plot(SVF,fraction=0.65,level=level)
title("zoomed out")
```

The first plot is zoomed in to the short lag behavior, while the second plot is zoomed out.
You can do this on the fly with `zoom(Cilla)` in R-studio. The variogram represents the average square distance traveled (vertical axis) within some time lag (horizontal axis).

For the long range behavior we can see that the variogram flattens (asymptotes) at approximately 20 days. This is, roughly, how coarse you need to make the timeseries so that methods assuming independence (no autocorrelation) can be valid. This includes, conventional kernel density estimation (KDE), minimum convex polygon (MCP), conventional species distribution modeling (SDM), and a host of other analyses.

The asymptote of our variogram is around 23 square km, and the fact that it takes roughly 20 days for the variogram to asymptote is indicative of the fact that the buffalo's location appears continuous at this timescale. This is also, roughly, the time it takes for the buffalo to cross its home range several times.

In the next sections we will "fit" models the variograms.
This is not a rigorous statistical fitting, but a good way of choosing candidate movement models and guessing at their parameters.
The variogram-fit parameter guestimates will then be fed into maximum likelihood estimation, which requires good initial guesses for non-linear optimization.

## Variogram Fitting the Hard Way

For pedagogical reasons we first "fit" models to the variogram by hand, while later we will use a much easier method.
we can manually guesstimate some continuous-time models for the aforementioned behavior with the commands

```{r,  fig.show='hold'}
m.iid <- ctmm(sigma=23 %#% "km^2")
m.ou <- ctmm(sigma=23 %#% "km^2",tau=6 %#% "day")
plot(SVF,CTMM=m.iid,fraction=0.65,level=level,col.CTMM="red")
title("Independent and identically distributed data")
plot(SVF,CTMM=m.ou,fraction=0.65,level=level,col.CTMM="purple")
title("Ornstein-Uhlenbeck movement")
```

where for both the models `m.iid` and `m.ou`, `sigma` is the asymptotic variance. In the Ornstein-Uhlenbeck (OU) model `m.ou`, `tau` is a single timescale that governs the autocorrelation in position and dictates the animal's home-range crossing time. The independent and identically distributed (IID) null model `m.iid` has no autocorrelation. Notice that all units are in meters and seconds, as output by the utility function `%#%` (see `help("%#%")`.

The IID model `m.iid` is obviously incorrect and in the zoomed in plot we can also see that the OU model `m.ou` is incorrectly linear at short lags, whereas the empirical variogram actually curves up for an hour or two before it becomes linear. Let us introduce a model that incorporates this behavior.

```{r,  fig.show='hold'}
m.ouf <- ctmm(sigma=23 %#% "km^2",tau=c(6 %#% "day",1 %#% "hour"))
plot(SVF,CTMM=m.ou,level=level,col.CTMM="purple",xlim=xlim)
title("Ornstein-Uhlenbeck movement")
plot(SVF,CTMM=m.ouf,level=level,col.CTMM="blue",xlim=xlim)
title("Ornstein-Uhlenbeck-F movement")
```

The confidence intervals at short lags are also very narrow, though both of these models look the same at coarser scales and so the discrepancy is only revealed by high resolution data.

```{r,  fig.show='hold'}
plot(SVF,CTMM=m.ou,fraction=0.65,level=level,col.CTMM="purple")
title("Ornstein-Uhlenbeck movement")
plot(SVF,CTMM=m.ouf,fraction=0.65,level=level,col.CTMM="blue")
title("Ornstein-Uhlenbeck-F movement")
```

The Ornstein-Uhlenbeck-F (OUF) model `m.ouf` introduces an additional autocorrelation timescale for the animal's velocity, so that it more closely matches the initial behavior of the variogram. The initial curve upwards tells us that there is continuity in the animal's velocity at this timescale. Conventional Markovian animal movement models do not capture this, which leads to the same kind of bias and underestimation of confidence intervals as when ignoring autocorrelation entirely.

The linear regime of the variogram (regular diffusion) is just as important as the asymptotic regime. In the linear regime it is reasonable to assume a Markovian model as with step selection functions (SSF) and Brownian bridges (BB). Therefore, the variogram has informed us as to how much we need to coarsen our data for it to be appropriate in many common analyses that neglect various aspects of movement.

## Variogram Fitting the Easy Way

The R-studio function `variogram.fit(SVF)` is much easier to use than guestimating the model parameters by hand as we did above. `variogram.fit` gives you sliders to choose the most visually appropriate parameters and save them to a global variable (`GUESS` by default).

### Variogram Error is Autocorrelated

It is important to note that variogram errors---the difference between the empirical variogram and the true semi-variance function---are themselves autocorrelated. Therefore, the smooth wiggling around that the point estimate does within the confidence bands is not necessarily meaningful. To demonstrate this, let us simulate some data from the OUF model `m.ouf` and look at its empirical variogram.

```{r,  fig.show='hold'}
# simulate fake buffalo with the same sampling schedule
willa <- simulate(m.ouf,t=Cilla$t)
plot(willa)
title("simulation")
# now calculate and plot its variogram
SVF2 <- variogram(willa)
plot(SVF2,CTMM=m.ouf,fraction=0.65,level=level,col.CTMM="blue")
title("simulation")
```

### Irregular Sampling Schedules

Random gaps in the data are acceptable and fully accounted for in both variogram estimation and model fitting.
However, if under any condition the sampling rate changes during data collection, then you will have to account for that in the variogram with the `dt` argument. In the following example, the collar was programmed to cycle between 1, 5, 25 hour sampling intervals

```{r, fig.show='hold'}
data("gazelle")
SVF3 <- variogram(gazelle[[18]])
plot(SVF3,fraction=0.85,level=level)
title("Default method")
# 1, 5, 25 hour sampling intervals
dt <- c(1,5,25) %#% "hour"
SVF3 <- variogram(gazelle[[18]],dt=dt)
plot(SVF3,fraction=0.85,level=level)
title("Multi method")
```

With small amounts of highly irregular data, you may also want to try `fast=FALSE`.

### Pooling Variograms

If multiple individuals exhibit similar movement behaviors, then we can pool their individual variograms to create a more precise population variogram. You should be careful though, if the individual movement behaviors and sampling schedules are not identical, then there will be discontinuities at lags where one timeseries runs out of data.

```{r, fig.show='hold'}
# buffalo 4 is bad
SVF4 <- lapply(buffalo[-4],variogram)
SVF4 <- mean(SVF4)
plot(SVF4,fraction=0.35,level=level)
title("Population variogram")
```

### Non-Stationarity in Stationary Variograms

Non-stationary behaviors, like a seasonal change in variance, is averaged over in the variogram. Moreover, if we fit a stationary model to non-stationary data, we are estimating an average effect. For instance, if an animal rests at night and diffuses at some rate `D` during the day, then without modeling the rest behavior we estimate an average of zero and `D`. Its not terribly detrimental to average over frequently repeated non-stationarity, but if an animal migrates once in a dataset then this behavior really needs to be in the model. These kinds of models will be included in future versions of `ctmm`.

Alternatively, you can also break up the data by hand into multiple behaviors and fit each behavior individually. This is highly recommended for migratory species, at present.

## More Models

### Telemetry Error

The GPS buffalo example do not exhibit telemetry errors that are significant enough to notice in our variograms. If we were working with ARGOS data or the high-resolution and unfiltered GPS data of a small animal, then we get a "nugget" effect that looks like an initial discontinuity at short time lags.

```{r,  fig.show='hold', echo=FALSE}
# ARGOS type errors
curve(1+x,0,5,xlab="Short time lag",ylab="Semi-variance",ylim=c(0,6))
points(c(0,0),c(0,1))
title("ARGOS")
# detector array type errors (qualitatively only)
curve((1-exp(-2*x))/(1-exp(-2/4)),0,1/4,xlab="Short time lag",ylab="Semi-variance",ylim=c(0,6),xlim=c(0,5),add=FALSE)
curve(3/4+x,1/4,5,xlab="Short time lag",ylab="Semi-variance",ylim=c(0,6),add=TRUE,xlim=c(0,5))
points(1/4,1)
title("Detector Array")
```

The height of this initial discontinuity corresponds to the variance of uncorrelated location errors. The second plot is the kind of initial discontinuity one has with detector array data. The end of the (slope) discontinuity is highlighted with a circle. This discontinuity is smooth because the movement and detection are correlated. The height of this initial discontinuity is also (at least roughly) the variance of the location errors.

If your tracking data is GPS and is annotated with a "DOP" column, then this information is automatically imported into the `telemetry` object and only has to be flagged in the `ctmm` model by setting `error=TRUE`. Because of some convenient mathematical relations, fitting without telemetry errors involves 1-4 fewer parameters and is, therefore, faster and more reliable than fitting with telemetry error. So if your data are not swamped by error and do not contain bursts of high-frequency fixes (which would also be relatively swamped by 5-10 meter GPS error), then you might want to apply a two stage fitting method, first from the variogram guess without telemetry error and then from that output with telemetry error. Here we demonstrate this method with Leroy the fisher cat from the `move` package:

```{r}
DATA <- as.telemetry(system.file("extdata","leroy.csv.gz",package="move"))
# default model guess
GUESS <- ctmm.guess(DATA,interactive=FALSE)
# first fit without telemetry error
FITS <- list()
FITS$NOERR <- ctmm.fit(DATA,GUESS)
# second fit based on first with telemetry error
GUESS <- FITS$NOERR
GUESS$error <- TRUE
FITS$ERROR <- ctmm.fit(DATA,GUESS)
# model improvement
summary(FITS)
```

Also, see `help(uere)` if you have the benefit of calibration data.

### Cycles and Periodicities

### One-Off Migrations

### Repeated Migrations

# Maximum Likelihood

## Maximum Likelihood Fitting the Hard Way 

Here we take our guestimates from variogram fitting the hard way and perform model selection manually.
Later we will finish everything off the easy way.
First let us fit each of our proposed models `m.iid`, `m.ou`, `m.ouf`, store the corresponding best-fit result in `M.IID`, `M.OU`, `M.OUF`, and then compare some of their outputs.

```{r,  fig.show='hold'}
M.IID <- ctmm.fit(Cilla,m.iid)
summary(M.IID)
```

```{r,  fig.show='hold'}
M.OU <- ctmm.fit(Cilla,m.ou)
summary(M.OU)
```

```{r,  fig.show='hold'}
M.OUF <- ctmm.fit(Cilla,m.ouf)
summary(M.OUF)
```

Notice how tiny the (Gaussian) area uncertainty is in IID model `M.IID`. Let us look into some details of the models.

```{r,  fig.show='hold'}
FITS <- list(IID=M.IID,OU=M.OU,OUF=M.OUF)
summary(FITS)
```

`AICc` is the (linearly) corrected Akaike information criteria. AIC balances likelihood against model complexity in a way that is good if we want to make optimal predictions. A lower AIC is better. Getting the AIC to go down by 5 is great, while getting the AIC to go down by 10 is awesome. Our AIC is going down by thousands.

The fit parameter `DOF[mean]` is the number of degrees of freedom worth of data we have to estimate the stationary mean parameter, assuming that the model is correct. Notice that the IID model perceives thousands of independent data points, while the autocorrelated OU and OUF models only see a handful of independent data points. This is why the IID model produced tiny confidence intervals on the predicted (Gaussian) area.

## Maximum Likelihood Fitting the Easy Way

If you have a complex, hypothetical model in mind, say the OUF `m.ouf` as you would get from variogram fitting the easy way, then you can perform model selection more conveniently with the `ctmm.select` function. `ctmm.select` considers the initial guess (hypothesis) and then iterates this model to select the best model based upon an information criteria.

```{r}
FITZ <- ctmm.select(Cilla,m.ouf,verbose=TRUE,level=1)
summary(FITZ)
```

Note that `level=1` is set here for demonstration purposes only, to attempt every nearest model no matter how unlikely it may seem.
The isotropic and anisotropic (`isotropic=FALSE`) flags correspond to circular and elliptical covariances respectively---an option we did not consider above.
<!-- The Ornstein-Uhlenbeck-F (OUF) model is continuous in position and velocity (`M.OUF`),
the Ornstein-Uhlenbeck (OU) model is continuous in position only (`M.OU`),
and the independent and identically distributed (IID) model has no autocorrelation (`M.IID`). -->
The IID model was never considered here by `ctmm.select` because it first requires selecting OU over OUF in the nested model hierarchy.
See `help("ctmm")` for more options.

## Variograms as a Disgnostic for Maximum Likelihood

Now its time to make sure that our selected model is explaining the most significant features of the animal's movement. Let us plot the variogram again with our fit models

```{r,  fig.show='hold'}
plot(SVF,CTMM=FITS,col.CTMM=c("red","purple","blue"),fraction=0.65,level=0.5)
title("zoomed out")
plot(SVF,CTMM=FITS,col.CTMM=c("red","purple","blue"),xlim=xlim,level=0.5)
title("zoomed in")
```

Notice that the purple OU model `M.OU` is significantly biased downward and is underestimating diffusion. This is because the continuous-velocity behavior at short time lags, which `M.OU` does not account for, is throwing off the estimate. The IID model `M.IID` is ignoring autocorrelation completely, while the OU model `M.OU` is ignoring autocorrelation in the buffalo's velocity.

While the OUF `M.OUF` is the selected model among all candidates, `M.OUF` looks slightly biased upwards in comparison to the variogram fit `m.ouf`. Some of this is due to sampling variability, which be (partially) remedied in the `variogram` method by increasing the `res` argument. Differences here can also arise from not accounting for telemetry error, which is unfortunately not annotated in this data.
