Introduction to argosTrack

Introduction

The argosTrack package is intended to model animal movement data collected by the Argos system with state-space models. It was originally constructed for the analyzes in Albertsen et al. (2015), which was based on the Continuous Time Correlated Random Walk (Johnson et al. 2008) movement model. Since then, the package has been extended to include several movement models, and is now based on reference classes. The new reference class based implementation gives a moduled based interface that is (hopefully) easier to use, maintain and extend. The package has three types of modules: Observation, Measurement, Movement and Animal. The Observation class is used to contain the data to be modelled. The Measurement class is used to specify the measurement part of the state-space model, which can be done (almost) independently of the actual observations. The Movement classes implements the process part of the state-space model, and finally the Animal class collects the parts to a full state-space model, which can be fitted.

Included data

The package includes two data sets from Albertsen et al. (2015); a track from a subadult ringed seal and a track from an adult ringed seal. The adult data is accessed by

library(argosTrack)
d0 <- adult_ringed_seal
head(d0)
##       id                date lc    lat     lon
## 1 106373 2011-10-30 13:33:20  B 56.405 -78.889
## 2 106373 2011-10-30 13:42:06  A 56.365 -78.847
## 3 106373 2011-10-30 14:18:19  0 56.396 -78.878
## 4 106373 2011-10-30 14:31:08  B 56.388 -78.863
## 5 106373 2011-10-30 15:04:19  1 56.400 -78.833
## 6 106373 2011-10-30 15:24:11  1 56.421 -78.847

Likewise, the subadult data is accessed by

d <- subadult_ringed_seal
head(d)
##      id                date lc    lat     lon
## 1 43847 2010-10-20 01:31:57  B 56.526 -79.261
## 2 43847 2010-10-20 02:13:06  1 56.583 -79.174
## 3 43847 2010-10-20 03:12:42  A 56.630 -79.178
## 4 43847 2010-10-20 03:54:43  1 56.638 -79.228
## 5 43847 2010-10-20 06:03:39  A 56.616 -79.126
## 6 43847 2010-10-20 07:23:00  B 56.646 -79.094

We will use the subadult data to illustrate the functionality of the package. For convenience, we remove the ‘Z’ location class, which is usually considered invalid.

d <- subset(d,lc != "Z")

Observation class

The first step in analysing the data using argosTrack is collecting the observation data in the Observation class. The input required is the longitudes, latitudes, time stamps and Argos location classes of the track data.

obs <- Observation(lon = d$lon,
                   lat = d$lat,
                   dates = as.POSIXct(d$date),
                   locationclass = d$lc
                   )
obs
## Observations:
## -------------
## 
## Number of observations: 3583 
## First date: 2010-10-20 
## Last date: 2011-04-06 
## Location classes:
## qual
##  GPS    3    2    1    0    A    B    Z 
##    0   92  209  392  263  854 1773    0

A plot of the track, latitudes over time, and longitudes over time can be produced by the plot function

plot(obs)
Default plot of movement track from subadult ringed seal.

Default plot of movement track from subadult ringed seal.

Each of the plot components can be produced by the functions plotMap, plotLat, and plotLon respectively. The plot functions can be customized by a list of base graphichs arguments as illustrated with the adult ringed seal data.

dateNum <- as.numeric(as.POSIXct(d0$date)) - as.numeric(as.POSIXct(d0$date[1]))
plotMap(Observation(lon = d0$lon,
                   lat = d0$lat,
                   dates = as.POSIXct(d0$date),
                   locationclass = d0$lc
                   ),args = list(col=as.numeric(d0$lc),
                                 cex = 3 * dateNum/tail(dateNum,1)
                                 )
        )
Map of movement track from adult ringed seal.

Map of movement track from adult ringed seal.

Measurement class

We now specify the measurement equation of the state-space model using the Measurement class. The main argument is model, which specifies whether the data should be considered normally distributed ("n") or t-distributed ("t"). In this example, we model the measurements by a normal distribution. Another important argument is nauticalObs, which specifies whether the observations should be transformed to nautical miles instead of latitudes and longitudes. The default (FALSE) is to model the latitudes and longitudes directly.

meas <- Measurement(model="n")
meas
## 
## 
## Measurement model
## -----------------
## 
## Measurement distribution: n 
## Use nautical observations: FALSE 
## Variance parameters:
##           GPS        3        2         1         0        A        B       Z
## Latitude    1 3.674437 108.4620 2214.5400 2113.2874 7211.267 69064.77 8886111
## Longitude   1 5.187231  60.7277  372.5607  442.5706 1620.030 61685.24 8886111

Movement class

The movement model is specified by Movement classes. The Movement class does not itself implement a movement model, and can not be used directly. Instead, each specific movement model in the package has its own class, which inherits the Movement class. Currently, the following classes implements movement models: CTCRW, DCRW, DIRAC, DSBHN, DSBW, GDCRW, OUL, RW. Information on each model can be seen in the documentation.

In this example, we use the RW class, which implements a random walk model. The class requires the time stamps at which the true position of the animal should be estimated; in our case, this is every (unique) time stamp of the data. By default the time steps are calculated in hours. This can be changed by the timeunit argument.

mov <- RW(unique(as.POSIXct(d$date)))
mov
## 
## 
## Movement model:
## -------------
## 
## Model: Random Walk (RW) 
## Movement parameters:
## numeric(0)
## Movement variance parameters:
## log(sigma_lat) log(sigma_lon) 
##              0              0 
## Number of latent variables: 6834 
## Using nautical states: FALSE 
## Using time unit: hours

Like the observations, the Movement classes can be plotted.

plot(mov)
Default plot of a non-fitted movement class.

Default plot of a non-fitted movement class.

Since we have not yet fitted the movement model to data, the plot currently shows the initial values.

Animal class

To collect the parts, the Animal class is used. The class requires an Observation, Movement, and Measurement model. Further, it allows us to give an identifier to the animal.

anim <- Animal(obs,mov,meas, as.character(d$id[1]))
anim
## 
## 
## Animal  43847 
## ------------------
## ------------------
## 
## 
## The model has not been fitted.
## 
## Observations:
## -------------
## 
## Number of observations: 3583 
## First date: 2010-10-20 
## Last date: 2011-04-06 
## Location classes:
## qual
##  GPS    3    2    1    0    A    B    Z 
##    0   92  209  392  263  854 1773    0 
## 
## 
## Measurement model
## -----------------
## 
## Measurement distribution: n 
## Use nautical observations: FALSE 
## Variance parameters:
##           GPS        3        2         1         0        A        B       Z
## Latitude    1 3.674437 108.4620 2214.5400 2113.2874 7211.267 69064.77 8886111
## Longitude   1 5.187231  60.7277  372.5607  442.5706 1620.030 61685.24 8886111
## 
## 
## Movement model:
## -------------
## 
## Model: Random Walk (RW) 
## Movement parameters:
## numeric(0)
## Movement variance parameters:
## log(sigma_lat) log(sigma_lon) 
##              0              0 
## Number of latent variables: 6834 
## Using nautical states: FALSE 
## Using time unit: hours

Using the plot function shows a combination of the data from the Observation class and the Movement class.

plot(anim)
Default plot of a non-fitted animal class.

Default plot of a non-fitted animal class.

Maximum likelihood estimation

Maximum likelihood estimation of the model using the Laplace approximation can be done with the fitTrack function; a wrapper for a combination of the TMB package and nlminb. The silent argument can be used to reduce the output written to the terminal during optimization. Depending on the movement model, further arguments to combine or fix parameters can be given. More information can be seen in the documentation for each movement model.

fitTrack(anim, silent=TRUE)

After estimating, the Animal object has been updated with the results.

anim
## 
## 
## Animal  43847 
## ------------------
## ------------------
## 
## 
## 
## Estimation results
## ------------------
## 
## Time to estimate: 11.573 
## Convergence message:
## relative convergence (4) 
## 
## Negative log likelihood: -8513.404 
## Maximum gradient component: 5.6e-04 
## Observations:
## -------------
## 
## Number of observations: 3583 
## First date: 2010-10-20 
## Last date: 2011-04-06 
## Location classes:
## qual
##  GPS    3    2    1    0    A    B    Z 
##    0   92  209  392  263  854 1773    0 
## 
## 
## Measurement model
## -----------------
## 
## Measurement distribution: n 
## Use nautical observations: FALSE 
## Variance parameters:
##                    GPS            3            2            1            0
## Latitude  1.679878e-06 6.172607e-06 9.429987e-13 7.759373e-06 0.0004674616
## Longitude 6.508404e-06 3.376059e-05 6.818810e-13 2.040193e-04 0.0024417375
##                      A           B        Z
## Latitude  0.0005639578 0.005589375 14.92758
## Longitude 0.0048599753 0.034029064 57.83440
## 
## 
## Movement model:
## -------------
## 
## Model: Random Walk (RW) 
## Movement parameters:
## numeric(0)
## Movement variance parameters:
## log(sigma_lat) log(sigma_lon) 
##      -3.923377      -3.239045 
## Number of latent variables: 6834 
## Using nautical states: FALSE 
## Using time unit: hours

This is also evident when plotting the object.

plot(anim)
Default plot of a fitted animal class.

Default plot of a fitted animal class.

Not only the Animal class, but also the Measurement and Movement class are updated.

mov
## 
## 
## Movement model:
## -------------
## 
## Model: Random Walk (RW) 
## Movement parameters:
## numeric(0)
## Movement variance parameters:
## log(sigma_lat) log(sigma_lon) 
##      -3.923377      -3.239045 
## Number of latent variables: 6834 
## Using nautical states: FALSE 
## Using time unit: hours

Because the random walk model does not have any movement parameters, this can not be seen when printing, but if we plot the object, the estimated locations are changed.

plot(mov,sd=FALSE)
Default plot of a fitted movement class.

Default plot of a fitted movement class.

The estimated track and its standard errors can be extracted with the function getTrack.

head(getTrack(anim))
##                 dates    id obs.lat obs.lon obs.lc  est.lat   est.lon
## 1 2010-10-20 01:31:57 43847  56.526 -79.261      B 56.58098 -79.17995
## 2 2010-10-20 02:13:06 43847  56.583 -79.174      1 56.58362 -79.17744
## 3 2010-10-20 03:12:42 43847  56.630 -79.178      A 56.61849 -79.19952
## 4 2010-10-20 03:54:43 43847  56.638 -79.228      1 56.63748 -79.21984
## 5 2010-10-20 06:03:39 43847  56.616 -79.126      A 56.63893 -79.15022
## 6 2010-10-20 07:23:00 43847  56.646 -79.094      B 56.66085 -79.11750
##        sd.lat       sd.lon
## 1   0.2260276    0.3002310
## 2 622.4857211 1336.7294643
## 3   0.2854371    0.3355757
## 4   0.2366838    0.3063652
## 5   0.2286369    0.3018959
## 6   0.2267473    0.3008016

Simulating data

To simulate from the state-space model (or any part of it), the simTrack function can be used. The following code simulates observations. The first argument is a Measurement class, specifying the model to simulate from; the second argument is the number of replications, while the third argument is an Observation class, which specifies the time points and location classes to use. The parameters used is the result from fitTrack. If fitTrack has not yet been used on the object, the initial parameters are used.

set.seed(123)
sobs <- simTrack(meas,1,obs)
plot(t(sobs[,,1]), pch = 16, col= rgb(0,0,0,0.3), xlab = "x", ylab = "y")
Plot of simulated observations

Plot of simulated observations

By plotting the simulated observations seperately for each location class, we clearly see how the variance increases with the location class.

layout(matrix(1:6,3))
sapply(2:7,function(i)
    plot(t(sobs[,as.numeric(obs$qual)==i,1]),
         pch = 16,
         xlab = "x", ylab = "y", main = levels(obs$qual)[i],
         col= rgb(0,0,0,0.3),
         ylim = range(sobs[2,,1]),
         xlim = range(sobs[1,,1]))
    )
Plot of simulated observations per location class

Plot of simulated observations per location class

Likewise, the function can be used on a Movement class. The first argument is the ‘Movement’ class to simulate from, the second argument is the number of replications, and the third argument is the initial position. If an initial position is not given, (0,0) is used.

set.seed(456)
smov <- simTrack(mov,1, c(0,0))
plot(t(smov[,,1]), xlab = "x", ylab = "y")
Plot of simulated movement

Plot of simulated movement

If the simTrack function is used on an Animal object, both movement and observations are simulated according to the full state-space model. Further, a new Animal object is returned. Note, that only the observations are changed from the object used in the simulation.

set.seed(789)
sa <- simTrack(anim,1)
## Save the new Animal object as an2
an2 <- sa[3,1]$Animal
an2
## 
## 
## Animal  43847_copy 
## ------------------
## ------------------
## 
## 
## 
## Estimation results
## ------------------
## 
## Time to estimate: 11.573 
## Convergence message:
## relative convergence (4) 
## 
## Negative log likelihood: -8513.404 
## Maximum gradient component: 5.6e-04 
## Observations:
## -------------
## 
## Number of observations: 3583 
## First date: 2010-10-20 
## Last date: 2011-04-06 
## Location classes:
## qual
##  GPS    3    2    1    0    A    B    Z 
##    0   92  209  392  263  854 1773    0 
## 
## 
## Measurement model
## -----------------
## 
## Measurement distribution: n 
## Use nautical observations: FALSE 
## Variance parameters:
##                    GPS            3            2            1            0
## Latitude  1.679878e-06 6.172607e-06 9.429987e-13 7.759373e-06 0.0004674616
## Longitude 6.508404e-06 3.376059e-05 6.818810e-13 2.040193e-04 0.0024417375
##                      A           B        Z
## Latitude  0.0005639578 0.005589375 14.92758
## Longitude 0.0048599753 0.034029064 57.83440
## 
## 
## Movement model:
## -------------
## 
## Model: Random Walk (RW) 
## Movement parameters:
## numeric(0)
## Movement variance parameters:
## log(sigma_lat) log(sigma_lon) 
##      -3.923377      -3.239045 
## Number of latent variables: 6834 
## Using nautical states: FALSE 
## Using time unit: hours
plot(an2)

The new Animal object can be used directly in fitTrack to update the parameter values and estimated locations.

fitTrack(an2, silent = TRUE)
## Warning in sqrt(diag(cov)): NaNs produced
## Warning in sqrt(diag(object$cov.fixed)): NaNs produced
## Warning in sqrt(as.numeric(object$diag.cov.random)): NaNs produced
## Warning in sqrt(diag(object$cov.fixed)): NaNs produced
## Warning in sqrt(as.numeric(object$diag.cov.random)): NaNs produced

Now, the new object is updated

an2
## 
## 
## Animal  43847_copy 
## ------------------
## ------------------
## 
## 
## 
## Estimation results
## ------------------
## 
## Time to estimate: 9.153 
## Convergence message:
## relative convergence (4) 
## 
## Negative log likelihood: -8630.438 
## Maximum gradient component: 7e-03 
## Observations:
## -------------
## 
## Number of observations: 3583 
## First date: 2010-10-20 
## Last date: 2011-04-06 
## Location classes:
## qual
##  GPS    3    2    1    0    A    B    Z 
##    0   92  209  392  263  854 1773    0 
## 
## 
## Measurement model
## -----------------
## 
## Measurement distribution: n 
## Use nautical observations: FALSE 
## Variance parameters:
##                    GPS            3            2            1            0
## Latitude  2.061595e-06 7.575201e-06 1.157266e-12 6.323384e-06 0.0004419754
## Longitude 3.233689e-06 1.677389e-05 3.387916e-13 1.975097e-04 0.0028722179
##                      A           B        Z
## Latitude  0.0005738976 0.005281115 18.31956
## Longitude 0.0046444013 0.033081094 28.73492
## 
## 
## Movement model:
## -------------
## 
## Model: Random Walk (RW) 
## Movement parameters:
## numeric(0)
## Movement variance parameters:
## log(sigma_lat) log(sigma_lon) 
##      -3.945934      -3.254290 
## Number of latent variables: 6834 
## Using nautical states: FALSE 
## Using time unit: hours
plot(an2)
Plot of estimated locations for the simulation

Plot of estimated locations for the simulation

The estimated parameters can easily be compared to the parameters used in the simulation. Below are the true and estimated variance parameters for the movement model.

anim$movement$varianceParameters
## [1] -3.923377 -3.239045
an2$movement$varianceParameters
## [1] -3.945934 -3.254290
## vcov is the covariance of both movement and variance parameters
2 * sqrt(diag(an2$movement$vcov))
## logSdState logSdState 
## 0.04896202 0.05564258

References

Albertsen, Christoffer Moesgaard, Kim Whoriskey, David Yurkowski, Anders Nielsen, and Joanna Mills Flemming. 2015. “Fast Fitting of Non-Gaussian State-Space Models to Animal Movement Data via Template Model Builder.” Ecology 96 (10): 2598–2604. https://doi.org/10.1890/14-2101.1.
Johnson, Devin S., Joshua M. London, Mary-Anne Lea, and John W. Durban. 2008. “Continuous-Time Correlated Random Walk Model for Animal Telemetry Data.” Ecology 89 (5): 1208–15. https://doi.org/10.1890/07-1032.1.