crwMLE: Fit Continuous-Time Correlated Random Walk Models to Animal Telemetry Data

Description

The function uses the Kalman filter to estimate movement paramters in a state-space version of the continuous-time movement model. Separate models are specified for movement portion and the location error portion. Each model can depend on time indexed covariates. A haul out model where movement is allowed to completely stop, as well as, a random drift model can be fit with this function.

Usage

crwMLE(mov.model=~1, err.model, stop.model=NULL, drift.model=FALSE,
       data, coord=c("x", "y"), polar.coord=TRUE, Time.name,
       initial.state, theta, fixPar, method="L-BFGS-B", control=NULL,
       initialSANN=NULL)

Arguments

mov.model

formula object specifying the time indexed covariates for movement parameters.

err.model

A 2-element list of formula objects specifying the time indexed covariates for location error parameters.

stop.model

formula object giving the covariate for the stopping portion of the model.

drift.model

logical indicating whether or not to include a random drift component.

data

data.frame object containg telemetry and covariate data.

coord

A 2-vector of character values giving the names of the "X" and "Y" coordinates in data.

polar.coord

logical indicating location are in latitude and longitude.

Time.name

character indicating name of the location time column

initial.state

list object containg the inital state of the Kalman filter.

theta

starting values for parameter optimization.

fixPar

Values of parameters which are held fixed to the given value.

method

Optimization method that is passed to optim.

control

Control list which is passed to optim.

initialSANN

Control list for optim when simulated annealing is used for obtaining start values. See details

Value

A list with the following elements:
parParameter maximum likelihood estimates
seStandard error of MLE
ci95% confidance intervals for parameters
CmatParameter covariance matrix
loglikMaximized log-likelihood value
aicModel AIC value
initial.stateIntial state provided to crwMLE for model fitting
coordCoordinate names provided for fitting
fixParFixed parameter values provided
convergenceIndicator of convergence (0 = converged)
messageMeesages given by optim during parameter optimization
stop.modelModel provided for stopping variable
random.driftLogical value indicating random drift model
mov.modelModel description for movement component
err.modelModel description for location error component
n.parnumber of parameters
nmsparameter names
y"latitude" locations
x"longitude" locations
n.movnumber of movement parameters
n.errXnumber or location error parameters for "longitude" error model
n.errYnumber or location error parameters for "latitude" error model
stop.mfcovariate for stop indication in stopping models
polar.coordCharacter vector giving coordinate names in data
runTimeTime used to fit model
initInitial values for parameter optimization

Details

A full model specification involves 4 components: a movement model, a stopping model, 2 location error models, and a drift indication. The movement model (mov.model) specifies how the movement parameters should vary over time. This is a function of specified, time-indexed, covariates. The movement parameters (sigma for velocity variation and beta for velocity autocorrelation) are both modeled with a log link as par = exp(eta), where eta is the linear predictor based on the covariates. The err.model specification is a list of 2 such models, one for longitude and one for latitude (in that order) location error. If only one location error model is given, it is used for both coordinates (parameter values as well). If drift.model is set to TRUE, then, 2 additional parameters are estimated for the drift process, a drift variance and a beta multiplier. The inital.state is a list with the following elemets (with the exact names): a1.y A vector with initial state values for the latitude coordinate. It have 2 elemets (location at time 1, velocity at time 1) for non-drift models and 3 elemets for drift models (location at time 1, velocity at time 1, drift velocity at time 1) for driftmodels, P1.y Covarince matrix for the state at time 1 (measure of uncertainty for your inital state) a1.y, a1.x Same as a1.y, but in the longitude coordinate, P1.x Same as P1.y, but in the longitude coordinate. theta and fixPar are vectors with the appropriate number or parameters. theta contains only those paraemters which are to be estimated, while fixPar contains all parameter values with NA for parameters which are to be estimated. The data set specified by data must contain a numeric column which is used as the time index for analysis. The column name is specified by the Time.name argument. Also, for stopping models, the stopping covariate must be between 0 and 1 inclusive, with 1 representing complete stop of the animal (no true movement, however, location error can still occur) and 0 represent unhindered movement. The coordinate location should have NA where no location is recorded, but there is a change in the movment covariates. The CTCRW models can be difficult to provide good initial values for optimization. If initialSANN is specified then simulated annealing is used first to obtain starting values for the specified optimaization method. If simulated annealing is used first, then the returned init list of the crwFit object will be a list with the results of the simulated annealing optimization.

Examples

Run this code

data(harborSeal)
head(harborSeal)
## Calculate Log multipliers for Argos error
argosClasses <- c("3", "2", "1", "0", "A", "B")
ArgosMultFactors <- data.frame(Argos_loc_class=argosClasses,
                               errX=log(c(1, 1.5, 4, 14, 5.21, 20.78)),
                               errY=log(c(1, 1.5, 4, 14, 11.08, 31.03)))
hsNew <- merge(harborSeal, ArgosMultFactors, by=c("Argos_loc_class"), all=TRUE)
hsNew <- hsNew[order(hsNew$Time), ]
head(hsNew)

## Initial state values
initial.dry <- list(
  a1.x=c(harborSeal$longitude[1],0),
  a1.y=c(harborSeal$latitude[1],0),
  P1.x=diag(c(1,1)),
  P1.y=diag(c(1,1))
)

##Fit model as given in Johnson et al. (2008) Ecology 89:1208-1215

fit1 <- crwMLE(
  mov.model=~1, err.model=list(x=~errX, y=~errY), stop.model=~DryTime,
  data=hsNew, coord=c("longitude","latitude"), Time.name="Time",
  initial.state=initial.dry, fixPar=c(NA, 1, NA, 1, NA, NA, NA),
  control=list(maxit=2000, trace=1, REPORT=1)
)

fit1
str(fit1)

##Use simulated annealing to obtain start value
fit2 <- crwMLE(
  mov.model=~1, err.model=list(x=~errX, y=~errY), stop.model=~DryTime,
  data=hsNew, coord=c("longitude","latitude"), Time.name="Time",
  initial.state=initial.dry, fixPar=c(NA, 1, NA, 1, NA, NA, NA),
  control=list(maxit=2000, trace=1, REPORT=1),
  initialSANN=list(maxit=100, temp=5, tmax=5, trace=1, REPORT=2)
)

##See simulated annealing start values
fit2$init