Joint (hierarchical) correlated random walk with a rotation-invariant process covariance

Source: R/jcrw_nll.R

Returns the RTMB objective function for a continuous-time correlated random walk fitted jointly to multiple individual tracks, with parameters shared among individuals according to a share_control specification.

jcrw_nll(dat)

Arguments

dat

a named list of model data assembled by jsfilter

Value

a function of the parameter list, suitable for RTMB::MakeADFun

Details

This is the R-language counterpart of src/TMB/sub/crw.hpp, extended across individuals and reparameterised so that the quantities shared among individuals do not depend on which way each animal travelled.

Why the reparameterisation. crw.hpp describes the velocity innovation covariance with D_x, D_y and rho_p, which are the axes and tilt of an ellipse measured against the projected coordinate axes. All three change when an animal travels a different bearing, even if its movement is otherwise identical. Pooling them across a deployment therefore averages quantities that mean different things for different animals, and an among-individual variance estimated from them is partly a statement about heading.

The parameterisation. The same covariance is written

  V = exp(m) * expm(A)        A = [ u1   u2 ]
                                  [ u2  -u1 ]

A is symmetric with zero trace, so exp(A) has determinant one and det(V) = exp(2m) exactly. A has eigenvalues +/- r, where r is sqrt(u1^2 + u2^2), so V has eigenvalues exp(m +/- r) and its major axis lies at theta = atan2(u2, u1) / 2. That gives

m

overall movement magnitude, half the log determinant. Rotation invariant: rotate the coordinate frame and it does not move.

r

log anisotropy, half the log ratio of the ellipse axes. Also rotation invariant.

theta

orientation. This is the only frame-dependent piece, and it is never pooled.

The 2 in atan2(u2, u1) / 2 is what makes u well behaved: an ellipse is unchanged by a half turn, so orientation lives on a half circle, and using 2 theta maps that onto a full circle. u is then an ordinary unconstrained vector with no wrapping to handle.

Because A is traceless and symmetric, exp(A) = cosh(r) I + sinh(r)/r A, so neither an eigendecomposition nor a matrix inverse is needed - the log determinant is 2m and the quadratic form is written out directly. That keeps the AD tape small, which matters more in RTMB than in C++.