---
title: "foieGras-basics"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{foieGras-basics}
  %\VignetteEngine{knitr::rmarkdown}
  \usepackage[utf8]{inputenc}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
library(foieGras)
library(dplyr)
library(ggplot2)
library(sf)
```

## Disclaimer
This vignette is an extended set of examples to highlight the `foieGras` package's functionality. Please, do NOT interpret these examples as instructions for conducting analysis of animal movement data. Numerous essential steps in a proper analysis have been left out of this document. It is your job to understand your data, ensure you are asking the right questions of your data, and that the analyses you undertake appropriately reflect those questions. We can not do this for you!

### foieGras models
This vignette provides a (very) brief overview of how to use `foieGras` to filter animal track locations obtained via the Argos satellite system or via processed light-level geolocation (GLS). `foieGras` provides two state-space models (SSM's) for filtering (ie. estimating "true" locations and associated movement model parameters, while accounting for error-prone observations):  

- a simple Random Walk model, `rw`
- a Correlated Random Walk model, `crw`  

Both models are continuous-time models, that is, they account for time intervals between successive observations, thereby naturally accounting for the irregularly-timed nature of most Argos data. We won't dwell on the details of the models here, those will come in a future paper, except to say there may be advantages to choosing one over the other in certain circumstances. The Random Walk model tends not to deal well with small to moderate gaps (relative to a specified time step) in observed locations and can over-fit to particularly noisy data. The Correlated Random Walk model can often deal better with these small to moderate data gaps and smooth through noisy data but tends to estimate nonsensical movement through larger data gaps.

Additionally, `foieGras` provides fast models (`mpm`, `jmpm`) for estimating a behavioural index along animals' tracks (see Jonsen et al. 2019 Ecology 100:e02566 for details). The `mpm` is fit to individual tracks, whereas the `jmpm` is fit to multiple individual track simultaneously with a variance parameter that is estimated jointly across the tracks. This latter model can often better resolve subtle changes in movement behaviour along tracks that lack much contrast in movements.

### input data
`foieGras` expects data to be provided in one of several possible formats.

1) a `data.frame` or `tibble` that looks like this
```{r data 1, echo = FALSE}
data(ellie, package = "foieGras")
head(ellie)
```
where the Argos data are provided via CLS Argos' Kalman filter model (KF) and include error ellipse information for each observed location.  

2) a `data.frame` or `tibble` that looks like this
```{r data 2, echo = FALSE}
data(ellies, package = "foieGras")
head(ellies)
```
where the Argos data are provided via CLS Argos' Least-Squares model (LS) and do not include error ellipse information.

3) a data.frame or tibble that includes observations with missing KF error ellipse information
```{r data 3, echo = FALSE}
data(ellie, package = "foieGras")
ellie[3:5, c("smaj","smin","eor")] <- NA
head(ellie)
```
in this situation, `foieGras` treats observations with missing error ellipse information as though they are LS-based observations.

4) an `sf` object where observations have any of the previous 3 structures and also include `CRS` information
```{r data 4, echo = FALSE, message=FALSE}
data(ellie, package = "foieGras")
foo <- sf::st_as_sf(ellie, coords=c("lon","lat"), crs = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs") 
foo <- sf::st_transform(foo, crs = "+proj=stere +lat_0=-90 +lat_ts=-71 +lon_0=0 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=km +no_defs")
head(foo)
```

5) a `data.frame`, `tibble` or `sf` object where processed GLS data are provided and include longitude and latitude error SD's (in degrees)
```{r data 5, echo = FALSE}

  tibble(
  id = rep(54632, 5),
  date = seq(Sys.time(), by = "12 hours", length.out = 5),
  lc = rep("G", 5),
  lon = seq(100, by = 0.5, length = 5),
  lat = seq(-55, by = 1, length = 5),
  lonerr = rexp(5, 1 / 0.5),
  laterr = rexp(5, 1 / 1.5)
  )
```

### fitting a foieGras model
model fitting for quality control of locations is comprised of 2 steps: a prefilter step where a number of checks are made on the input data (see `?foieGras::prefilter` for details), including applying the `argsofilter::sdafilter` to identify extreme outlier observations. Additionally, if the input data are not supplied as an `sf` object, `prefilter` guesses at an appropriate projection (typically world mercator, EPSG 3395) to apply to the data. The SSM is then fit to this projected version of the data. Users invoke this process via the `fit_ssm` function:
```{r fit_ssm, message=FALSE}
## load foieGras example data
data(ellie)
## prefilter and fit Random Walk SSM using a 24 h time step
fit <- fit_ssm(ellie, model = "rw", time.step = 24, verbose = 0)
```
these are the minimum arguments required: the input data, the model ("rw" or "crw") and the time.step (in h) to which locations are predicted. Additional control can be exerted over the prefiltering step, via the `vmax`, `ang`, `distlim`, `spdf` and `min.dt` arguments. see `?foieGras::fit_ssm` for details, the defaults for these arguments are quite conservative, usually leading to relative few observations being flagged to be ignored by the SSM. Additional control over the SSM fitting step can also be exerted but these should rarely need to be accessed by users and will not be dealt with here.

### accessing and visualizing model fit objects
Simple summary information about the `foieGras` fit can be obtained by calling the fit object:
```{r fit summary}
fit$ssm[[1]]
```
and a summary `plot` method allows a quick visual of the SSM fit to the data:
```{r fit plot, fig.width=7,fig.height=5}
# plot time-series of the predicted values
plot(fit, what = "predicted", type = 1)
plot(fit, what = "fitted", type = 2)
```
The predicted values (red) are the state estimates predicted at regular time intervals, specified by `time.step` (here every 24 h). These estimates are plotted on top of the observations that passed the prefilter stage (blue points and blue rug at bottom). Fitted values are the state estimates corresponding to the time of each observation; their time series are plotted by default - `plot(fit)`. A 2-D time series plot of the track is invoked by the argument `type = 2`.

As SSMs are latent variable models, evaluating their goodness of fit is less straightforward than simpler statistical models without latent variables. We can use One-Step-Ahead (prediction) residuals via `foieGras::osar`. Here we use `osar` to compare SSM fits of the `rw` and `crw` model to the same example southern elephant seal data.
```{r osar comp, fig.width=7, fig.height=5, message=FALSE}

## fit crw SSM
fitc <- fit_ssm(ellie, model = "crw", time.step = 24, verbose = 0)

## calculate OSA resids for both models
fit_res <- osar(fit)
fitc_res <- osar(fitc)

## plot residuals
plot(fit_res)
plot(fitc_res)
```
The `crw` model appears provide a better fit than the `rw` model, with standardized OSA residuals conforming more closely to a theoretical Normal distribution. One note of caution when calculating OSA residuals, the underlying `TMB::oneStepPredict` method is currently experimental and can require *considerable* computation time, especially when calculating across multiple individual fits.  

Estimated tracks can be mapped using the `fmap` function, which uses the `foieGras`-applied projection (Global Mercator). Projections can be changed easily via the `crs` argument in the form of a proj4string (as in the example, below).
```{r ggplot map, fig.width=5, fig.height=5, message=FALSE}
## map ssm-predicted values without observations
fmap(fitc, what = "predicted", obs = FALSE)

## change projection to Antarctic Polar Stereographic centred on 
##  the approximate mid-point of the track
fmap(fitc, what = "predicted", 
     crs = "+proj=stere +lat_0=-90 +lat_ts=-71 +lon_0=85 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=km +no_defs")
```

The estimated locations can be accessed for further analysis, custom mapping, etc... by using the `grab` function. They can be returned as a projected `sf` tibble or as a simple unprojected tibble. Note, that for all `foieGras` outputs the `x`, `y`, `x.se` and `y.se` units are in km.
```{r grab}
## grab fitted locations from fit object as a projected sf object 
plocs_sf <- grab(fitc, what = "f")

## grab predicted locations in unprojected form, returning as a tibble
plocs <- grab(fitc, "p", as_sf = FALSE)

## unprojected form looks like this
plocs
```

`fit_ssm` can be applied to single tracks as shown, it can also fit to multiple individual tracks in a single input `tibble` opr `data.frame`. The SSM is fit to each individual separately. The resulting output is a compound `tibble` with rows corresponding to each individual `foieGras` fit object. The `converged` column indicates whether each model fit converged successfully. 
```{r multi-fits}
# load 2 southern elephant seal example data
data(ellies)

fit2 <- fit_ssm(ellies, vmax = 10, model = "crw", time.step = 48, verbose = 0)

# list fit outcomes for both seals
fit2
```
individual `id` is displayed in the 1st column, all fit output (`ssm`) in the 2nd column, and `convergence` status of each model fit is displayed in the 3rd column

`fmap` automatically handles ssm fit objects with multiple individuals, plotting all on a single map
```{r ggplot map 2, fig.width=5, fig.height=5, message=FALSE}
## map predicted values and observations
fmap(fit2, "p", obs = TRUE)
```

A behavioural index can be estimated from locations provided they occur regularly in time and they either have minimal location error (i.e. GPS data) or they have been ssm filtered. We can fit the `mpm` to `foieGras` ssm-predicted locations. Here we use the ssm fits to the s. elephant seal data
```{r fit mpm, message=FALSE}
## fit mpm separately to each individual track
fmp <- fit2 %>% 
  grab(., "p", as_sf = FALSE) %>%
  select(id, date, lon, lat) %>%
  fit_mpm(., model = "mpm")

fmp
```

We can visualize the estimated behavioural index (move persistence) as a time series for each seal. The move persistence parameter $g_t$ ranges continuously from 0 (little persistence, indicative of area-restricted movements) to 1 (high persistence, indicative of directed movements).
```{r plot mpm ts 1, fig.width=7, fig.height=5, message=FALSE}
## plot mpm estimates by individual seal
grab(fmp, "fitted") %>% 
  ggplot() +
  geom_point(aes(date, g, colour = g)) +
  scale_colour_viridis_c(limits = c(0,1)) +
  ylim(0,1) +
  facet_wrap(~ id, scales = "free_x", ncol = 1)

```

A joint move persistence model `jmpm` is also available for fitting to multiple individuals. This model fits to all individuals simultaneously, estimating a joint random walk variance parameter that can often better resolve subtle variations in $g_t$. 

We can explore the spatio-temporal variation in movement behaviour by plotting the $g_t$ values along each seal's track, but first we have to merge the ssm-predicted locations with the move persistence estimates using `foieGras::join()`
```{r plot mpm track, fig.width=7, fig.height=5, message=FALSE}

## join ssm predicted locations and move persistence values together
fmp_locs <- join(fit2, fmp, as_sf = FALSE)

ggplot(fmp_locs) +
  geom_point(aes(lon, lat, colour = g)) +
  scale_colour_viridis_c(limits = c(0,1))
```


