Simulates multiple individual tracks from the jcrw process
and observation models with known parameter values, for testing whether
fit_ssm recovers them.
sim_jcrw(
A = 20,
n = 200,
dt = 6,
irregular = TRUE,
D = 0.005,
sd_m = 0.5,
aniso = 2,
orient = NULL,
error = c("LS", "KF"),
tau = c(1, 1),
rho_o = 0,
psi = 1,
lc.prob = c(`3` = 0.05, `2` = 0.1, `1` = 0.15, `0` = 0.2, A = 0.2, B = 0.3),
crs = paste("+proj=laea +lat_0=-57 +lon_0=70", "+datum=WGS84 +units=km +no_defs"),
start = c(0, 0),
start.date = "2020-01-01 00:00:00"
)number of individuals
number of locations per individual; a scalar, or a vector of
length A
mean time interval between locations, in hours
logical; draw intervals from a gamma distribution with
mean dt (TRUE, the default) rather than using a regular interval
population mean diffusion magnitude, on the same scale as crw's
D_x and D_y. The default gives plausible elephant seal speeds at the
default track length; see Details
standard deviation of log(2 * D) among individuals. Zero
gives identical individuals
ratio of the axes of the velocity innovation covariance; a
scalar shared by all individuals, or a vector of length A. 1 is isotropic
orientation of the major axis in degrees, a vector of length
A. NULL (default) draws orientations uniformly
observation error model: "LS" (default) for Argos
least-squares classes, "KF" for error ellipses
scaling on the error multiplication factors, length 2
correlation between the x and y observation errors, "LS"
only
scaling on the ellipse semi-minor axis, "KF" only
named vector of Argos location class probabilities
projection to simulate and return in
two-element vector giving the starting coordinate, in the
units of crs
start date-time, as a character string or POSIXt
an sf object in the requested projection with columns id,
date, lc and, for error = "KF", smaj, smin and eor. The
generating parameters and the latent states are attached as the attribute
"truth", for comparison against a fit
Tracks are simulated directly in a projected coordinate frame and
returned as an sf object carrying that projection. fit_ssm respects an
sf input, so nothing is re-projected between simulating and fitting and
no map distortion enters the comparison. That matters here: on a Mercator
grid the apparent diffusion magnitude is inflated by sec(latitude)^2, so
simulating in longitude and latitude would confound the recovery test with
the projection.
The process is the one in jcrw_nll:
with each individual's velocity innovation covariance built from its magnitude, anisotropy and orientation:
where m = log(2 * D), u1 = rho * cos(2 * theta),
u2 = rho * sin(2 * theta) and rho = log(aniso) / 2. V then has
eigenvalues exp(m +/- rho), so the ratio of its axes is aniso and its
major axis lies at theta.
Individual magnitudes are drawn as m[i] ~ N(log(2 * D), sd_m^2). Set
sd_m well away from zero to test whether the hierarchical model can
recover among-individual variation - a fit to data with no such variation
can only ever show sd_m collapsing to its lower bound, which tests
nothing.
A note on choosing D. Velocity in a correlated random walk is an
integrated random walk, so it has no stationary distribution: its standard
deviation grows as sqrt(2 * D * T) over a track of total duration T.
Simulated speeds therefore increase with track length, and a D taken from
a real fit will not generally produce plausible speeds when run forward -
a real fit is held in check by the observations, forward simulation is not.
At six-hourly intervals over 300 steps, D = 0.005 gives a median speed of
about 2.4 km/h, in the right range for a southern elephant seal, while
D = 0.05 gives about 7 km/h and D = 0.5 around 25 km/h. Pick D for
the speeds you want at the track length you are simulating, and check with
sqrt(colSums(attr(d, "truth")$states[[1]]$v^2)). This is a property of
the process model, not of the simulator.
Argos error ellipse parameters are drawn from the same log-normal and von
Mises distributions, fitted to real Argos data, that sim uses.
Note that the errors here are generated with eor converted to radians;
the internal ellp.par applies sin() and cos() to eor in degrees,
so its errors do not match the ellipse parameters recorded beside them.
## 20 individuals differing markedly in movement magnitude
d <- sim_jcrw(A = 20, n = 200, D = 0.005, sd_m = 0.6, aniso = 2)
truth <- attr(d, "truth")
round(truth$D, 3)
#> [1] 0.009 0.002 0.006 0.007 0.010 0.008 0.004 0.001 0.003 0.009 0.003 0.002
#> [13] 0.015 0.017 0.005 0.004 0.004 0.008 0.005 0.007