library(ggplot2)
library(patchwork)
library(nlme)
te_paper <- "#f5f4ee"
te_ink <- "#16241d"
te_body <- "#2c3a31"
te_forest <- "#275139"
te_rust <- "#b5534e"
te_gold <- "#c9b458"
te_line <- "#dad9ca"
theme_datasheet <- function() {
theme_minimal(base_size = 12) +
theme(plot.background = element_rect(fill = te_paper, colour = NA),
panel.background = element_rect(fill = te_paper, colour = NA),
panel.grid.major = element_line(colour = te_line, linewidth = 0.3),
panel.grid.minor = element_blank(),
text = element_text(colour = te_body),
plot.title = element_text(colour = te_ink, face = "bold"),
plot.subtitle = element_text(colour = te_body),
axis.text = element_text(colour = te_body),
strip.text = element_text(colour = te_ink))
}
class_cols <- c(resident = te_forest, nomad = te_line, dispersal = te_gold, migrant = te_rust)Classifying migrants from net squared displacement
R
movement ecology
migration
model selection
simulation
ecology tutorial
Choosing net squared displacement curves by AIC labels stationary home-range residents as migrants more often the slower they roam. Floors and AR(1) fits in R.
Corrected 14 October 2026: the curve search, Nelder-Mead from three grid shapes and a single line search for the resident rate, often stopped short of the best least squares fits, and the nls fits that beat it break the migrant’s timing constraint; the search now starts from many more shapes, seeds the disperser and migrant fits from each other and fits the resident rate on a fine grid, and every share, figure and comparison has been recomputed with it. The resident shares fall or stay the same under every rule, and with autocorrelated errors and AIC no resident is kept at fifteen days, where the post had said that this fit kept as many as independent errors and moved the mislabelled residents from the migrant to the disperser label; it now labels them mostly as migrants on short trips, and at fifteen days the choice is worse than with independent errors.
Corrected 8 October 2026: the post described the resident curve fit as a Nelder-Mead polish from grid starts and said that the nls check sized how far the search stopped short of the best migrant fit, so that it only slightly understated the false migrant rates; it now says that the resident curve is fitted by a line search and that the check does not size the shortfall of the disperser and migrant fits, which, if anything, makes the resident shares too high.
A herd of red deer winters low in a mountain valley, and twenty hinds carry GPS collars fitted in January. Each collar takes one fix a day for a year. Some of the hinds will climb to summer range in May and come back in the autumn, some will stay put, and the report needs a number: the proportion of the population that migrates. The standard way to get it is to turn each animal-year into a single series, the net squared displacement from the first fix, fit a small family of curves to it and let an information criterion pick the curve. A plateau is a resident, a single step is a disperser, a step up and a step back down is a migrant, a straight line is a nomad.
That model set is Bunnefeld and colleagues’ 2011 approach, and the migrateR package of Spitz, Hebblewhite and Stephenson packages it with a set of constraints on the fitted parameters. Singh, Allen and Ericsson wrote a paper of clarifications and caveats for the approach, including that short-distance migrants can be classified as residents or nomads, and Cagnacci and colleagues compared it with other ways of assessing migration on the same data. None of that is new here. What this post measures is the other side of the error: how often a hind that never left her range is called a migrant, and what sets the rate.
This site has three posts that touch the pieces. Correlated random walks and net displacement uses net squared displacement as a diffusion null, and because a single walk is noisy it averages the curve over a thousand replicate walks before reading it. A collar gives one walk per animal per year, and the classification is made on that one noisy series. Home ranges in R: MCP versus kernel density simulates a range as an Ornstein-Uhlenbeck process and shows that the autocorrelation time, roughly how long the animal takes to forget where it was, sets how much a track can say about the range. That time scale turns out to be the whole story below. Brownian motion or Ornstein-Uhlenbeck? counts how often an information criterion prefers the richer model on data from the simpler one, on phylogenies; the same kind of count is made here on tracks, where the richer models are the migrant and disperser curves.
The residents are simulated as a stationary Ornstein-Uhlenbeck range, so the true answer is always resident and every other label is an error. The post fits the four curves carefully, counts how the choice depends on the range’s autocorrelation time, and then tries three repairs: a stricter criterion, the displacement floor that migrateR offers, and autocorrelated errors in the curve fit.
Four curves for one animal-year
day <- 0:364 # one daily fix per day, first fix on day 0
n_fix <- length(day)
s_km <- 1 # stationary SD of the range on each axis, km
tau_set <- c(1, 3, 7, 15, 30, 60) # position autocorrelation time, days
n_res <- 200 # resident animal-years per tau
tau_mig <- 15 # tau for the migrant runs
dist_set <- c(1, 2, 4, 8) # migration distance, km
n_mig <- 200 # migrant animal-years per distance
leave_day <- 120; return_day <- 270; move_phi <- 5
floor_set <- c(3, 5) # displacement floors, km
tau_ar <- c(3, 15, 60) # tau for the AR(1) refits
n_ar <- 200 # residents refitted per tau
r95_km <- sqrt(qchisq(0.95, 2)) * s_km # radius holding 95 per cent of positions
logistic <- function(theta, phi) 1 / (1 + exp((theta - day) / phi))
# 2-D OU around a fixed centre by its exact daily update, plus an optional
# seasonal shift of the centre along x; NSD is measured from the first fix
sim_nsd <- function(tau, s = s_km, dist = 0) {
a <- exp(-1 / tau)
ou_axis <- function() as.numeric(stats::filter(rnorm(n_fix, 0, s * sqrt(1 - a^2)), a,
method = "recursive", init = rnorm(1, 0, s)))
x <- ou_axis() + dist * (logistic(leave_day, move_phi) - logistic(return_day, move_phi))
y <- ou_axis()
(x - x[1])^2 + (y - y[1])^2
}A resident hind is a two-dimensional Ornstein-Uhlenbeck process around a fixed centre: on each axis her position is pulled back towards the centre with time constant tau, and her stationary distribution has a standard deviation of 1 km on each axis, so 95 per cent of her positions lie within 2.45 km of the centre. Fleming and colleagues treat this time scale as the range crossing time, roughly how long the animal takes to cross its range. It runs from 1 to 60 days in the design, from an animal that crosses its range every day or two to one that spends weeks in one corner before drifting to another. A migrant is the same process whose centre moves 1 to 8 km along one axis in a logistic step around day 120 and back around day 270. The daily update is exact, so the tracks carry no discretisation error, and the first fix is drawn from the stationary distribution, so the collar starts wherever the animal happens to be. All of these constants, and the numbers of simulated animal-years, were fixed before the first run.
The four curves are the ones in Bunnefeld and colleagues’ paper, written for net squared displacement (NSD) against day. The nomad is a straight line through the origin. The resident rises to a plateau delta as one minus a negative exponential. The disperser is a logistic step of height delta at day theta with width phi. The migrant is a logistic step up at theta minus a step of the same height down at theta2, each with its own width, and, as in the migrateR parameterisation, the midpoint of the return lies at least twice the sum of the two widths after the midpoint of the departure. The mixed migrant curve, which returns to a different place, is left out.
Fitting the curves without failed fits
Fitting these curves with nls needs starting values, and a fit that stops with an error quietly becomes a missing model, which changes the choice. The fits here avoid that by the shape of the models. Every curve is delta times a shape that does not involve delta, so for any fixed shape the least squares delta has a closed form, truncated at zero. The fit is then a search over the shape parameters alone. The one-parameter resident curve is fitted by a fine grid over its rate and a line search around the best grid value. The disperser and migrant curves can have several local optima, so Nelder-Mead runs from many starting shapes and is restarted from its own result until it stops improving: the best shapes of a starting grid, the best of a wider grid with steps before, inside and after the collar year (among them returns in the last days of the year and trips of a few days), the best sharp steps between two days, and each other’s fit, because a migrant whose return lies far beyond the year is a disperser. The lowest residual sum of squares wins. Nothing in that search can stop with an error, so no animal-year is ever missing a model.
n_par <- c(nomad = 1, resident = 2, dispersal = 3, migrant = 5) # excluding sigma
# curve shapes without delta; p on the optimiser scale (log rates and log phi)
shape_of <- function(model, p) switch(model,
resident = 1 - exp(-exp(p[1]) * day),
dispersal = logistic(p[1], exp(p[2])),
migrant = logistic(p[1], exp(p[3])) - logistic(p[2], exp(p[4])))
valid_par <- function(model, p) model != "migrant" ||
p[2] >= p[1] + 2 * exp(p[3]) + 2 * exp(p[4])
rho_set <- c(0, seq(0.05, 0.995, length.out = 60))
# least squares for delta >= 0 given the shape; with ar = TRUE the errors are
# AR(1), fitted by conditional least squares with rho profiled over rho_set
profile_delta <- function(X, y, ar = FALSE) {
X <- as.matrix(X)
if (!ar) {
xy <- crossprod(X, y)[, 1]; xx <- colSums(X^2); dl <- pmax(xy / xx, 0)
return(list(rss = sum(y^2) - 2 * dl * xy + dl^2 * xx, delta = dl, rho = rep(0, ncol(X))))
}
X1 <- X[-1, , drop = FALSE]; X0 <- X[-n_fix, , drop = FALSE]
y1 <- y[-1]; y0 <- y[-n_fix]
c11 <- crossprod(X1, y1)[, 1]; c10 <- crossprod(X1, y0)[, 1]
c01 <- crossprod(X0, y1)[, 1]; c00 <- crossprod(X0, y0)[, 1]
s11 <- colSums(X1^2); s10 <- colSums(X1 * X0); s00 <- colSums(X0^2)
q11 <- sum(y1^2); q10 <- sum(y1 * y0); q00 <- sum(y0^2)
rss <- rep(Inf, ncol(X)); dl <- rho <- rep(0, ncol(X))
for (r in rho_set) {
xy <- c11 - r * (c10 + c01) + r^2 * c00
xx <- s11 - 2 * r * s10 + r^2 * s00
dd <- pmax(xy / xx, 0)
rr <- q11 - 2 * r * q10 + r^2 * q00 - 2 * dd * xy + dd^2 * xx
better <- rr < rss
rss[better] <- rr[better]; dl[better] <- dd[better]; rho[better] <- r
}
list(rss = rss, delta = dl, rho = rho)
}
# starting grids on the optimiser scale, one shape per column
start_grid <- list(
resident = cbind(log(exp(seq(log(0.003), log(3), length.out = 40)))),
dispersal = as.matrix(expand.grid(seq(5, 360, 5), log(c(1, 3, 7, 15, 30)))),
migrant = local({
g <- as.matrix(expand.grid(seq(10, 350, 10), seq(10, 350, 10), log(c(2, 6, 15)), log(c(2, 6, 15))))
g[apply(g, 1, function(p) valid_par("migrant", p)), ]
}))
start_shape <- lapply(names(start_grid), function(m)
apply(start_grid[[m]], 1, function(p) shape_of(m, p)))
names(start_shape) <- names(start_grid)
# one shape at a time, with the sums over the track computed once per track
rss_shape <- function(x, trk, ar) {
if (!ar) {
xy <- sum(x * trk$y); xx <- sum(x * x); dl <- max(xy / xx, 0)
return(c(trk$yy - 2 * dl * xy + dl^2 * xx, dl, 0))
}
x1 <- x[-1]; x0 <- x[-n_fix]; r <- rho_set
xy <- sum(x1 * trk$y1) - r * (sum(x1 * trk$y0) + sum(x0 * trk$y1)) + r^2 * sum(x0 * trk$y0)
xx <- sum(x1 * x1) - 2 * r * sum(x1 * x0) + r^2 * sum(x0 * x0)
dl <- pmax(xy / xx, 0); rr <- trk$q - 2 * dl * xy + dl^2 * xx; i <- which.min(rr)
c(rr[i], dl[i], r[i])
}
# wider starting sets for the search; the grid-only comparison below keeps start_grid,
# and the search draws no random numbers
wide_grid <- list(
resident = cbind(seq(log(1e-4), log(20), length.out = 4000)),
dispersal = as.matrix(expand.grid(seq(-200, 900, 10), log(c(0.5, 1, 2, 4, 8, 15, 30, 60, 120, 250)))),
migrant = local({
g <- rbind(as.matrix(expand.grid(seq(-50, 800, 25), seq(0, 1000, 50), log(c(1, 4, 15, 50, 150)), log(c(1, 4, 15, 50, 150)))),
as.matrix(expand.grid(-990, seq(-300, 450, 25), log(1), log(c(1, 4, 15, 50, 150, 400)))))
g[apply(g, 1, function(p) valid_par("migrant", p)), ]
}))
# migrant families the main grid covers poorly: a return within a few days of the end
# of the collar year, a slow decline from a return far beyond it, and short trips
late_grid <- local({
g <- as.matrix(expand.grid(seq(-250, 350, 25), c(361.5, 362.5, 363.5, 364.5, 366, 370),
log(c(0.25, 1, 4, 15, 50, 150, 300)), log(c(0.25, 1, 3))))
g[apply(g, 1, function(p) valid_par("migrant", p)), c(1, 2, 3, 4)]
})
slow_grid <- local({
g <- as.matrix(expand.grid(seq(-50, 350, 10), c(600, 800, 1000), log(c(0.25, 1, 4, 15)), log(c(150, 250, 400))))
g[apply(g, 1, function(p) valid_par("migrant", p)), ]
})
hump_grid <- local({
g <- as.matrix(expand.grid(seq(0, 360, 2), c(2, 3, 5, 8, 12, 20, 35, 60), log(c(0.5, 1.5))))
g <- cbind(g[, 1], g[, 1] + g[, 2], g[, 3], g[, 3])
g[apply(g, 1, function(p) valid_par("migrant", p)), ]
})
wide_grid$late <- late_grid; wide_grid$slow <- slow_grid; wide_grid$hump <- hump_grid
wide_shape <- lapply(names(wide_grid), function(m)
apply(wide_grid[[m]], 1, function(p) shape_of(if (m %in% c("late", "slow", "hump")) "migrant" else m, p)))
names(wide_shape) <- names(wide_grid)
for (m in names(wide_grid)) { # drop shapes that are zero on every day of the year
keep <- colSums(wide_shape[[m]]^2) > 0
wide_grid[[m]] <- wide_grid[[m]][keep, , drop = FALSE]
wide_shape[[m]] <- wide_shape[[m]][, keep, drop = FALSE]
}
# Nelder-Mead, restarted from its own result until it stops improving; for the
# migrant it alternates with a run in coordinates in which every point keeps the timing
# constraint p2 >= p1 + 2 phi1 + 2 phi2 (p2 = p1 + 2 phi1 + 2 phi2 + s^2)
to_q <- function(p) c(p[1], sqrt(max(p[2] - p[1] - 2 * exp(p[3]) - 2 * exp(p[4]), 0)), p[3], p[4])
to_p <- function(q) c(q[1], q[1] + 2 * exp(q[3]) + 2 * exp(q[4]) + q[2]^2, q[3], q[4])
polish <- function(par, objective, migrant = FALSE, max_runs = 30) {
val <- objective(par)
obj_q <- function(q) objective(to_p(q))
for (run in seq_len(max_runs)) {
old <- val
o <- optim(par, objective, control = list(maxit = 400))
if (o$value < val) { par <- o$par; val <- o$value }
if (migrant) {
o <- optim(to_q(par), obj_q, control = list(maxit = 400))
if (o$value < val) { par <- to_p(o$par); val <- o$value }
}
# a run on the offsets from the current point: optim's first simplex is then a
# tenth of a unit in every direction instead of a tenth of the largest parameter
p0 <- par
o <- optim(0 * p0, function(d) objective(p0 + d), control = list(maxit = 400))
if (o$value < val) { par <- p0 + o$par; val <- o$value }
if (run >= 2 && old - val <= 1e-10 * val) break
}
par
}
# the best shape in each block of timings, so that the starts are spread out
spread_best <- function(rss, grid, width, k) {
cell <- if (ncol(grid) <= 2) floor(grid[, 1] / width) else paste(floor(grid[, 1] / width), floor(grid[, 2] / width))
o <- order(rss); first <- o[!duplicated(cell[o])]
first[seq_len(min(k, length(first)))]
}
# sharp steps between two days: the best windows found by cumulative sums, as starts
# with a narrow width (the fit is almost flat in the timing of a sharp step); with AR(1)
# errors the windows are scored on the series y[t] - rho y[t - 1] for every rho of the grid
step_starts <- function(model, y, k, ar = FALSE) {
n <- length(y)
if (model == "dispersal") {
a <- 1:(n - 1) # on from day a
} else {
ab <- which(upper.tri(matrix(0, n, n)), arr.ind = TRUE)
a <- ab[, 1] - 1; b <- ab[, 2] - 1 # away on days a to b - 1
}
if (!ar) {
cs <- c(0, cumsum(y)); yy <- sum(y^2)
S <- if (model == "dispersal") cs[n + 1] - cs[a + 1] else cs[b + 1] - cs[a + 1]
len <- if (model == "dispersal") n - a else b - a
r <- yy - pmax(S, 0)^2 / len
w <- log(0.25)
} else {
r <- Inf
for (rho in rho_set) {
z <- y[-1] - rho * y[-n]; cz <- c(0, cumsum(z)); Q <- sum(z^2) # z[t] belongs to day t
if (model == "dispersal") {
xy <- z[a] + (1 - rho) * (cz[n] - cz[a + 1]); xx <- 1 + (1 - rho)^2 * (n - 1 - a)
} else {
za <- ifelse(a == 0, 0, z[pmax(a, 1)])
xy <- za + (1 - rho) * (cz[b] - cz[pmax(a, 1) + (a > 0)]) - rho * z[b]
xx <- (a > 0) + (1 - rho)^2 * (b - 1 - pmax(a, 1) + (a == 0)) + rho^2
}
r <- pmin(r, Q - pmax(xy, 0)^2 / xx)
}
w <- log(0.05)
}
# the best three also with a width of one day, and with the first (last) day of the
# window taking a part of the step (timings a quarter day either side, width 0.15)
part <- log(0.15)
if (model == "dispersal") {
i <- order(r)[seq_len(k)]
return(c(lapply(a[i], function(aa) c(aa - 0.5, w)), lapply(a[i[1:3]], function(aa) c(aa - 0.5, 0)),
unlist(lapply(a[i[1:3]], function(aa) lapply(c(-0.25, 0, 0.25), function(cc) c(aa + cc, part))), recursive = FALSE)))
}
cell <- paste(a %/% 20, b %/% 20)
o <- order(r); o <- o[!duplicated(cell[o])][seq_len(k)]
shifts <- list(c(-0.25, -0.25), c(0, 0), c(0.25, 0.25), c(-0.25, 0.25), c(0.25, -0.25))
c(lapply(o, function(j) c(if (a[j] == 0) -5 else a[j] - 0.5, b[j] - 0.5, w, w)),
lapply(o[1:3], function(j) c(if (a[j] == 0) -5 else a[j] - 0.5, b[j] - 0.5, 0, 0)),
unlist(lapply(o[1:3], function(j) lapply(shifts, function(cc)
c(if (a[j] == 0) -5 else a[j] + cc[1], b[j] - 1 + cc[2] + 0.5, part, part))), recursive = FALSE))
}
# resident: a fine grid over the log rate, then a line search between the neighbours
# of the best grid value; sigmoid curves: Nelder-Mead from the three best shapes of
# the starting grid, the spread-out best of the wider grid (and of the migrant
# families), the best sharp steps and any seeds, keeping the lowest residual sum of squares
fit_curve <- function(model, y, ar = FALSE, n_start = 3, n_wide = 12, n_step = 6, seeds = list(), grids = TRUE) {
trk <- list(y = y, yy = sum(y^2), y1 = y[-1], y0 = y[-n_fix])
trk$q <- sum(trk$y1^2) - 2 * rho_set * sum(trk$y1 * trk$y0) + rho_set^2 * sum(trk$y0^2)
if (model == "nomad") {
v <- rss_shape(day, trk, ar)
return(list(rss = v[1], delta = v[2], rho = v[3], par = numeric(0), top = v[2] * max(day)))
}
objective <- function(p) {
if (!valid_par(model, p) || any(abs(p) > 1000)) return(.Machine$double.xmax)
rss_shape(shape_of(model, p), trk, ar)[1]
}
if (model == "resident") {
g <- wide_grid$resident[, 1]
i <- which.min(profile_delta(wide_shape$resident, y, ar)$rss)
o <- optimize(objective, g[c(max(i - 1, 1), min(i + 1, length(g)))])
starts <- list(g[i], o$minimum)
} else {
starts <- seeds
if (grids) {
pr <- profile_delta(start_shape[[model]], y, ar)
pw <- profile_delta(wide_shape[[model]], y, ar)
iw <- spread_best(pw$rss, wide_grid[[model]], if (model == "migrant") 100 else 50, n_wide)
starts <- c(starts,
lapply(order(pr$rss)[seq_len(n_start)], function(i) start_grid[[model]][i, ]),
lapply(iw, function(i) wide_grid[[model]][i, ]),
step_starts(model, y, n_step, ar))
if (model == "migrant") for (fam in c("late", "slow", "hump")) {
pf <- profile_delta(wide_shape[[fam]], y, ar)
jf <- spread_best(pf$rss, wide_grid[[fam]][, 1, drop = FALSE], 100, 3)
starts <- c(starts, lapply(jf, function(i) wide_grid[[fam]][i, ]))
}
}
starts <- Filter(function(p) objective(p) < .Machine$double.xmax, starts)
if (length(starts) == 0) return(NULL)
starts <- lapply(starts, polish, objective = objective, migrant = model == "migrant")
}
vals <- sapply(starts, objective)
best_par <- starts[[which.min(vals)]]
best <- rss_shape(shape_of(model, best_par), trk, ar)
list(rss = best[1], delta = best[2], rho = best[3], par = best_par,
top = best[2] * max(shape_of(model, best_par)))
}
# with its return at the edge of the parameter box the migrant curve is close to the
# disperser (for a step narrow enough), and the disperser is the migrant's first step:
# each fit seeds the other until the disperser stops improving
fit_all <- function(y, ar = FALSE) {
fits <- lapply(c(nomad = "nomad", resident = "resident"), fit_curve, y = y, ar = ar)
d <- fit_curve("dispersal", y, ar)
m <- fit_curve("migrant", y, ar)
for (round in 1:5) {
m_new <- fit_curve("migrant", y, ar, seeds = list(c(d$par[1], 1000, d$par[2], 0)), grids = FALSE)
if (!is.null(m_new) && m_new$rss < m$rss) m <- m_new
d_new <- fit_curve("dispersal", y, ar, seeds = list(m$par[c(1, 3)]), grids = FALSE)
if (!is.null(d_new) && d_new$rss < d$rss * (1 - 1e-12)) d <- d_new else break
}
c(fits, list(dispersal = d, migrant = m))
}
# one column of results per animal-year; trip is the time between the migrant's two steps
year_row <- function(fits, ar) {
fits <- fits[names(n_par)]
c(sapply(fits, `[[`, "rss"), sapply(fits, `[[`, "delta"), sapply(fits, `[[`, "top"), ar = ar,
trip = unname(fits$migrant$par[2] - fits$migrant$par[1]))
}
fit_year <- function(y, ar = FALSE) year_row(fit_all(y, ar), ar)
# the choice; floor_km removes dispersal and migrant fits with delta below
# floor_km^2 (migrateR's mdelta); nomad_rule also removes the nomad fit when its
# largest fitted NSD is within floor_km^2 of the resident's (migrateR's mnr);
# parsimony takes the fewest parameters within 2 AIC units (migrateR's a.rule)
choose_model <- function(fr, pen = 2, floor_km = 0, nomad_rule = FALSE, parsimony = FALSE) {
ar <- fr[["ar"]]; m_obs <- n_fix - ar
ic <- m_obs * log(fr[1:4] / m_obs) + pen * (n_par + 1 + ar)
keep <- c(TRUE, TRUE, fr[7:8] >= floor_km^2)
if (nomad_rule) keep[1] <- abs(fr[[10]] - fr[[9]]) >= floor_km^2
ic[!keep] <- Inf
pick <- which.min(ic)
if (parsimony) pick <- which(ic - min(ic) <= 2)[1] # n_par is in increasing order
names(n_par)[pick]
}The same device gives a fit with AR(1) errors. For a fixed shape and a fixed autocorrelation rho the conditional least squares problem is again linear in delta, so rho is profiled over a grid of values from zero to 0.995 inside the same search, and the criterion counts rho as one more parameter.
The choice itself is a function of the fitted residual sums of squares and the two migrateR rules. Its displacement floor drops a disperser or migrant fit whose delta is below a minimum; since delta is in squared km, a floor of 3 km means delta below 9. Nomads have no delta, so the floor never touches them; migrateR has a separate rule for them, which drops the nomad fit when its largest fitted NSD differs from the largest fitted NSD of the resident curve by less than a set amount, here the same squared floor. By default migrateR also takes the model with the fewest parameters among those within 2 AIC units of the best, which is included as a third variant of the criterion.
set.seed(4417)
chk_y <- replicate(40, sim_nsd(30))
chk_nls <- t(apply(chk_y, 2, function(y) {
out <- c(); fa <- fit_all(y)
for (m in c("dispersal", "migrant")) {
pr <- profile_delta(start_shape[[m]], y); i <- which.min(pr$rss); p0 <- start_grid[[m]][i, ]
dat <- data.frame(dday = day, nsd = y)
fm <- if (m == "dispersal") {
try(nls(nsd ~ delta / (1 + exp((theta - dday) / phi)), dat,
start = list(delta = pr$delta[i], theta = p0[[1]], phi = exp(p0[[2]]))), silent = TRUE)
} else {
try(nls(nsd ~ delta / (1 + exp((theta - dday) / phi)) - delta / (1 + exp((theta2 - dday) / phi2)), dat,
start = list(delta = pr$delta[i], theta = p0[[1]], theta2 = p0[[2]],
phi = exp(p0[[3]]), phi2 = exp(p0[[4]]))), silent = TRUE)
}
nls_rss <- if (inherits(fm, "try-error")) NA else deviance(fm)
out <- c(out, nls_rss, fa[[m]]$rss)
}
# does the migrant nls fit keep the timing constraint of the migrant curve?
cf <- if (inherits(fm, "try-error")) NULL else coef(fm)
c(out, if (is.null(cf)) NA else cf[["phi"]] > 0 && cf[["phi2"]] > 0 &&
cf[["theta2"]] >= cf[["theta"]] + 2 * cf[["phi"]] + 2 * cf[["phi2"]])
}))
nls_fail <- colSums(is.na(chk_nls[, c(1, 3)]))
nm_gain <- n_fix * log(chk_nls[, c(1, 3)] / chk_nls[, c(2, 4)]) # AIC units, NA if nls failed
nls_beats <- colSums(nm_gain < -0.01, na.rm = TRUE) # nls found a lower RSS
nm_worst <- -apply(nm_gain, 2, min, na.rm = TRUE) # largest shortfall, AIC units
nm_ahead <- apply(nm_gain, 2, max, na.rm = TRUE) # largest gain over nls, AIC units
nls_out <- sum(nm_gain[, 2] < -0.01 & chk_nls[, 5] == 0, na.rm = TRUE) # lower, but outside the model
n_chk <- ncol(chk_y)
# the AR(1) profile against nlme::gnls with corAR1, dispersal curve, ten tracks
gn_out <- t(sapply(1:10, function(j) {
pf <- fit_all(chk_y[, j], ar = TRUE)$dispersal
fm <- try(gnls(nsd ~ delta / (1 + exp((theta - dday) / phi)),
data.frame(dday = day, nsd = chk_y[, j]),
start = c(delta = pf$delta, theta = pf$par[[1]], phi = exp(pf$par[[2]])),
correlation = corAR1(form = ~ dday)), silent = TRUE)
if (inherits(fm, "try-error")) return(c(NA, pf$rho, NA, pf$delta))
c(coef(fm$modelStruct$corStruct, unconstrained = FALSE)[[1]], pf$rho, coef(fm)[["delta"]], pf$delta)
}))
gn_fail <- sum(is.na(gn_out[, 1]))
gn_rho_gap <- max(abs(gn_out[, 1] - gn_out[, 2]), na.rm = TRUE)
gn_delta_gap <- max(abs(gn_out[, 3] / gn_out[, 4] - 1), na.rm = TRUE)
# scale: the same random numbers with a range half and twice as wide
set.seed(918); y_half <- replicate(20, sim_nsd(15, s = 0.5))
set.seed(918); y_twice <- replicate(20, sim_nsd(15, s = 2))
fr_half <- apply(y_half, 2, fit_year); fr_twice <- apply(y_twice, 2, fit_year)
same_class <- mean(apply(fr_half, 2, choose_model) == apply(fr_twice, 2, choose_model))
delta_ratio <- range(fr_twice[8, ] / fr_half[8, ])The checks come before the results. The first compares the search with nls, started from the best grid shape, on 40 resident tracks with tau of 30 days. For the disperser nls failed on 2 tracks and found a lower residual sum of squares than the search on 0 of the others (counting only gaps above a hundredth of an AIC unit), while the search came out lower by up to 38.4 AIC units. For the migrant nls failed on 21 of the 40 tracks. Where it converged it found a lower residual sum of squares on 6 tracks, by at most 2.6 AIC units, but on 6 of them its fit breaks the timing constraint of the migrant curve, which nls does not enforce and the search keeps; elsewhere the search came out lower by up to 34.8 AIC units.
The second fits the AR(1) disperser with nlme::gnls and a corAR1 structure on the first ten of the same tracks, starting from the profile fit. gnls stopped with an error on 8 of them. On the rest its rho differed from the profile by at most 0.011 and its delta by at most 9 per cent, a gap plausibly due to the conditional least squares used here against the full likelihood in gnls.
The last check is on scale rather than on the optimiser. Multiplying every position by a constant c multiplies every NSD by c squared and every residual sum of squares by c to the fourth power, which leaves every AIC difference unchanged. The same random numbers with a range half as wide and twice as wide gave the same class for 100 per cent of 20 tracks, and delta ratios of 16 to 16. Without a floor the choice depends on tau and on the migration distance over the range size, never on the range size alone. With a floor it depends on the floor over the range size, which matters for the floor section below.
Residents turn into migrants as tau grows
set.seed(20260917)
res_y <- lapply(tau_set, function(tau) replicate(n_res, sim_nsd(tau)))
res_fr <- lapply(res_y, function(ym) apply(ym, 2, fit_year))
share_tab <- function(fr_list, keys, rule, ...) {
do.call(rbind, lapply(seq_along(fr_list), function(i) {
cl <- apply(fr_list[[i]], 2, choose_model, ...)
data.frame(key = keys[i], class = names(n_par), rule = rule,
share = as.numeric(table(factor(cl, levels = names(n_par)))) / ncol(fr_list[[i]]))
}))
}
get_share <- function(tab, rule_name, cls) tab$share[tab$rule == rule_name & tab$class == cls]
mc_se <- function(p, nn) sqrt(p * (1 - p) / nn)
res_share <- rbind(share_tab(res_fr, tau_set, "AIC"),
share_tab(res_fr, tau_set, "AIC, parsimony rule", parsimony = TRUE),
share_tab(res_fr, tau_set, "BIC", pen = log(n_fix)))
aic_ok <- get_share(res_share, "AIC", "resident")
par_ok <- get_share(res_share, "AIC, parsimony rule", "resident")
bic_ok <- get_share(res_share, "BIC", "resident")
aic_mig <- get_share(res_share, "AIC", "migrant")
aic_dis <- get_share(res_share, "AIC", "dispersal")
nomad_max <- max(res_share$share[res_share$class == "nomad"])
se_max <- max(mc_se(c(aic_ok, par_ok, bic_ok), n_res))
# the same tracks classified from the starting grid alone, without the search
grid_year <- function(y) {
fits <- lapply(names(n_par), function(m) {
pr <- if (m == "nomad") profile_delta(day, y) else profile_delta(start_shape[[m]], y)
i <- which.min(pr$rss); c(pr$rss[i], pr$delta[i])
})
c(sapply(fits, `[`, 1), sapply(fits, `[`, 2), rep(0, 4), ar = 0)
}
i15 <- match(15, tau_set)
grid_ok <- mean(apply(apply(res_y[[i15]], 2, grid_year), 2, choose_model) == "resident")With AIC, the share of residents classified as residents was 0.750 at tau of 1 day, 0.430 at 3 days, 0.100 at 15 days and 0.020 at 60 days, from 200 animal-years per value with Monte Carlo standard errors of 0.035 or less. At tau of 60 days the migrant curve won for 0.715 of the residents and the disperser for 0.265. Nomads were almost never chosen: the largest nomad share in any cell, under any of the three criteria, was 0.005.
The error at the short end is not zero either. At a tau of one day, positions a day apart keep a correlation of only 0.37, and still 0.250 of those residents were given a sigmoid curve. A plausible reason is that the timing parameters of the sigmoid curves act like change points, which can buy more fit than the 2 AIC units each is charged, and that NSD measured from one fix is a skewed series whose spread changes with its level, for which a least squares criterion is only an approximation.
The parsimony rule of migrateR raised the resident share to 0.155 at 15 days, and BIC, whose penalty per parameter is 5.9 for a year of daily fixes, to 0.305. BIC kept 0.985 of the residents at 1 day and 0.840 at 3 days, but by 60 days it keeps only 0.055 of the residents. No penalty that is fixed per parameter can follow tau, because the extra fit a wandering resident offers the migrant curve grows with tau while the penalty does not.
The optimiser matters in the direction nobody would pick. Classifying the same 200 residents at tau of 15 days from the starting grid alone, without the search, kept 0.345 of them as residents against 0.100 with it. A coarse fit to the flexible curves flatters the method; the better the curves are fitted, the more false migrants appear.
plot_share <- res_share
plot_share$class <- factor(plot_share$class, levels = c("migrant", "dispersal", "nomad", "resident"))
plot_share$tau <- factor(plot_share$key, levels = tau_set)
ggplot(plot_share, aes(tau, share, fill = class)) +
geom_col(width = 0.75, colour = te_paper, linewidth = 0.2) +
facet_wrap(~ rule) +
scale_fill_manual(values = class_cols, name = "label") +
scale_y_continuous(breaks = seq(0, 1, 0.25)) +
labs(x = "range autocorrelation time tau (days)", y = "share of residents",
title = "Stationary residents labelled by curve choice",
subtitle = sprintf("%d animal-years per bar, daily fixes for a year", n_res)) +
theme_datasheet() +
theme(legend.position = "bottom")
example_panel <- function(y, lab) {
fa <- fit_all(y); pick <- choose_model(year_row(fa, FALSE))
fc <- fa[[pick]]
shape <- if (pick == "nomad") day else shape_of(pick, fc$par)
list(pts = data.frame(day = day, nsd = y, panel = lab),
fit = data.frame(day = day, nsd = fc$delta * shape, panel = lab, pick = pick))
}
ex <- list(example_panel(res_y[[1]][, 1], sprintf("resident, tau %d d", tau_set[1])),
example_panel(res_y[[5]][, 1], sprintf("resident, tau %d d", tau_set[5])))
set.seed(33)
ex[[3]] <- example_panel(sim_nsd(tau_mig, dist = 4), sprintf("migrant 4 km, tau %d d", tau_mig))
ex_pts <- do.call(rbind, lapply(ex, `[[`, "pts"))
ex_fit <- do.call(rbind, lapply(ex, `[[`, "fit"))
ex_lab <- unique(ex_fit[, c("panel", "pick")])
ex_pts$panel <- factor(ex_pts$panel, levels = unique(ex_pts$panel))
ex_fit$panel <- factor(ex_fit$panel, levels = levels(ex_pts$panel))
ex_lab$panel <- factor(ex_lab$panel, levels = levels(ex_pts$panel))
ggplot(ex_pts, aes(day, nsd)) +
geom_point(size = 0.6, colour = te_body, alpha = 0.45) +
geom_line(data = ex_fit, aes(colour = pick), linewidth = 1.1) +
geom_label(data = ex_lab, aes(x = 5, y = Inf, label = paste("AIC choice:", pick)),
hjust = 0, vjust = 1.3, size = 3.4, fill = te_paper, colour = te_ink) +
facet_wrap(~ panel, scales = "free_y") +
scale_colour_manual(values = class_cols, guide = "none") +
labs(x = "day of the collar year", y = "NSD from the first fix (km squared)",
title = "One animal-year, one curve") +
theme_datasheet()
The example panels show why the flexible curves win. A resident with a long autocorrelation time spends weeks in one part of her range before drifting to another, so her distance from the first fix wanders in slow swings instead of scattering around a plateau. In the middle panel it climbs through most of the year and drops near the end, and a wide step up with a late step down describes that better than any plateau.
A displacement floor, and what it costs
A floor asks for a minimum size of the step before a sigmoid curve may win. The scale check above says what a floor in km really is: a floor in units of the range. For this range, 95 per cent of positions lie within 2.45 km of the centre, so the floors of 3 and 5 km are 1.2 and 2.0 times that radius; the same floors on a range twice as wide would be half as strict. The migrants were run at tau of 15 days, 200 animal-years per distance.
floor_res <- rbind(
share_tab(res_fr, tau_set, "no floor"),
share_tab(res_fr, tau_set, "3 km floor", floor_km = floor_set[1]),
share_tab(res_fr, tau_set, "5 km floor", floor_km = floor_set[2]),
share_tab(res_fr, tau_set, "5 km floor and nomad rule", floor_km = floor_set[2], nomad_rule = TRUE))
f3_ok <- get_share(floor_res, "3 km floor", "resident")
f5_ok <- get_share(floor_res, "5 km floor", "resident")
f5_nom <- get_share(floor_res, "5 km floor", "nomad")
f5n_ok <- get_share(floor_res, "5 km floor and nomad rule", "resident")
set.seed(170926)
mig_y <- lapply(dist_set, function(dd) replicate(n_mig, sim_nsd(tau_mig, dist = dd)))
mig_fr <- lapply(mig_y, function(ym) apply(ym, 2, fit_year))
floor_mig <- rbind(
share_tab(mig_fr, dist_set, "no floor"),
share_tab(mig_fr, dist_set, "3 km floor", floor_km = floor_set[1]),
share_tab(mig_fr, dist_set, "5 km floor", floor_km = floor_set[2]),
share_tab(mig_fr, dist_set, "5 km floor and nomad rule", floor_km = floor_set[2], nomad_rule = TRUE))
m0_hit <- get_share(floor_mig, "no floor", "migrant")
m3_hit <- get_share(floor_mig, "3 km floor", "migrant")
m5_hit <- get_share(floor_mig, "5 km floor", "migrant")
m5_res <- get_share(floor_mig, "5 km floor", "resident")
m5n_hit <- get_share(floor_mig, "5 km floor and nomad rule", "migrant")Among residents at tau of 15 days, where AIC alone kept 0.100, the 3 km floor kept 0.735 and the 5 km floor 0.930. At 60 days the numbers were 0.020, 0.560 and 0.755. The floor removes sigmoid fits, and some of those animal-years then fall to the nomad line instead of the resident curve: with the 5 km floor the nomad share at 60 days was 0.110, against almost none without it. Adding the nomad rule at the same distance put the resident share back to 0.865. A floor that is not also applied to the nomad model moves part of the error to a different wrong label.
The cost falls on short migrants. Without a floor AIC labelled 0.820 of the 2 km migrants as migrants, 0.985 of the 4 km migrants and 1.000 of the 8 km migrants. The 3 km floor left 0.450, 0.950 and 1.000; the 5 km floor left 0.050, 0.495 and 1.000, and the lost 4 km migrants went to resident (0.500). A fitted delta is a squared distance measured from a first fix that is itself somewhere in the winter range, so it scatters around the square of the true distance, and a floor near the true distance removes a large part of the migrants at that distance.
The 1 km migrants are a warning about the whole table rather than about the floor. Their step is well inside the winter range, and without a floor 0.740 of them were called migrants, against 0.705 of the stationary residents at the same tau. At that distance the migrant label says almost nothing about whether the animal moved.
rule_cols <- c("no floor" = te_rust, "3 km floor" = te_gold,
"5 km floor" = te_forest, "5 km floor and nomad rule" = te_ink)
res_line <- floor_res[floor_res$class == "resident", ]
mig_line <- floor_mig[floor_mig$class == "migrant", ]
p_res <- ggplot(res_line, aes(key, share, colour = rule)) +
geom_line(linewidth = 0.8) + geom_point(size = 1.8) +
scale_x_log10(breaks = tau_set) +
scale_colour_manual(values = rule_cols, breaks = names(rule_cols), name = NULL) +
coord_cartesian(ylim = c(0, 1)) +
labs(x = "tau (days, log scale)", y = "residents labelled resident", title = "Residents kept") +
theme_datasheet()
p_mig <- ggplot(mig_line, aes(key, share, colour = rule)) +
geom_line(linewidth = 0.8) + geom_point(size = 1.8) +
scale_x_log10(breaks = dist_set) +
scale_colour_manual(values = rule_cols, breaks = names(rule_cols), name = NULL) +
coord_cartesian(ylim = c(0, 1)) +
labs(x = "migration distance (km, log scale)", y = "migrants labelled migrant",
title = "Migrants found") +
theme_datasheet()
(p_res + p_mig) + plot_layout(guides = "collect") +
plot_annotation(theme = theme_datasheet()) & theme(legend.position = "bottom")
Autocorrelated errors do not restore the choice
The statistical complaint about the curve fits is that daily NSD values are treated as independent. The obvious repair keeps the curves and gives the errors an AR(1) structure, as gnls with corAR1 would. The residents at three values of tau were refitted that way, 200 animal-years each, so the comparison with the independent-error fits is paired.
ar_idx <- match(tau_ar, tau_set)
ar_fr <- lapply(ar_idx, function(i) apply(res_y[[i]][, seq_len(n_ar)], 2, fit_year, ar = TRUE))
iid_fr <- lapply(ar_idx, function(i) res_fr[[i]][, seq_len(n_ar)])
ar_share <- rbind(share_tab(iid_fr, tau_ar, "independent errors, AIC"),
share_tab(iid_fr, tau_ar, "independent errors, BIC", pen = log(n_fix)),
share_tab(ar_fr, tau_ar, "AR(1) errors, AIC"),
share_tab(ar_fr, tau_ar, "AR(1) errors, BIC", pen = log(n_fix - 1)))
iid_ok <- get_share(ar_share, "independent errors, AIC", "resident")
ar_ok <- get_share(ar_share, "AR(1) errors, AIC", "resident")
arb_ok <- get_share(ar_share, "AR(1) errors, BIC", "resident")
iidb_ok <- get_share(ar_share, "independent errors, BIC", "resident")
ar_mig <- get_share(ar_share, "AR(1) errors, AIC", "migrant")
ar_dis <- get_share(ar_share, "AR(1) errors, AIC", "dispersal")
arb_nom <- get_share(ar_share, "AR(1) errors, BIC", "nomad")
ar_se <- max(mc_se(c(iid_ok, iidb_ok, ar_ok, arb_ok), n_ar))
# median time between the two steps of the migrant curves that AIC chose, at the middle tau
trip_of <- function(fr, ...) median(fr["trip", apply(fr, 2, choose_model, ...) == "migrant"])
ar_trip <- trip_of(ar_fr[[2]]); iid_trip <- trip_of(iid_fr[[2]])At tau of 3 days the AR(1) fit kept 0.495 of the residents against 0.430 for the same tracks with independent errors, and with BIC 0.935 (independent errors with BIC: 0.840). At 15 days it kept 0.000 against 0.100; the migrant share rose to 0.925 and the disperser share fell to 0.075. At 60 days the resident share was 0.005 against 0.020. BIC with AR(1) errors kept 0.145 and 0.025 at the two longer time scales, against 0.305 and 0.055 with independent errors and BIC, so the AR(1) structure makes the choice worse at 15 days and does not improve it at 60, and it sent 0.130 and 0.135 of the residents to the nomad line. Monte Carlo standard errors are 0.035 or less.
The AR(1) structure helps at the short time scale (with BIC, from 0.840 to 0.935) and not at the long ones, where its resident shares are lower than with independent errors under both criteria. Most of the gain of AR(1) errors with BIC over independent errors with AIC comes from the penalty, not from the error structure. When rho is close to one, the conditional least squares criterion compares mostly day-to-day changes, and the migrant curves it picks are mostly short trips rather than seasonal ones: at 15 days the migrant curves chosen by AIC had a median of 18 days between their two steps with AR(1) errors, against 279 days with independent errors. NSD from an Ornstein-Uhlenbeck range is also not an AR(1) series with constant variance; it is a squared distance whose spread grows with its level. A model that describes the positions themselves, a stationary range against a shifting one in continuous time, is the natural next step and is not tried here.
ar_plot <- ar_share
ar_plot$class <- factor(ar_plot$class, levels = c("migrant", "dispersal", "nomad", "resident"))
ar_plot$tau <- factor(sprintf("tau %d d", ar_plot$key), levels = sprintf("tau %d d", tau_ar))
ar_plot$rule <- factor(ar_plot$rule, levels = unique(ar_share$rule))
ggplot(ar_plot, aes(rule, share, fill = class)) +
geom_col(width = 0.7, colour = te_paper, linewidth = 0.2) +
facet_wrap(~ tau) +
scale_fill_manual(values = class_cols, name = "label") +
scale_x_discrete(labels = c("iid\nAIC", "iid\nBIC", "AR(1)\nAIC", "AR(1)\nBIC")) +
labs(x = NULL, y = "share of residents",
title = "AR(1) errors help only at the shortest tau",
subtitle = sprintf("%d paired animal-years per tau", n_ar)) +
theme_datasheet() +
theme(legend.position = "bottom")
What to report
Report the curve set, the criterion and every constraint as numbers: the floor on delta in km and the range size it was compared against, whether the nomad model was held to a matching rule, and whether the parsimony rule was on. A floor stated in km means nothing without the range size, because without a floor the choice depends on the range size only through tau and the migration distance over the range size.
Estimate the autocorrelation time of the resident ranges from the positions before the classification, and say it. A share of migrants from a population whose residents cross their range in a day is not comparable with the same share from a population whose residents take a month, even with identical collars and identical code.
Simulate the population’s own residents. A stationary Ornstein-Uhlenbeck range with the estimated range size and tau, run through the same fitting and choice code, gives the false migrant rate for that study directly; the code above is enough for that. If the rate is not small, report the migrant share next to it, or treat the classification as a screen for animals whose tracks are then read by eye.
Keep a short-distance migrant category open. Where the migration distance is not well beyond the range radius, none of the rules measured here separates them from residents at tau of 15 days or more without losing most of them, and a floor large enough to protect the residents removes those migrants by construction.
Honest limits
Every resident here has a single stationary range with an isotropic, Gaussian Ornstein-Uhlenbeck process, and every migrant makes one symmetric trip with fixed dates. Real winter ranges drift, are shaped by terrain and snow, and have more than one centre; each of those adds structure that the sigmoid curves can take for movement, so there is no reason to expect real residents with the same tau to be misclassified less often than these.
The fixes are daily and perfect: no location error, no gaps, no four-hourly schedule. More fixes per day at the same tau add points without adding independent information, so a fixed penalty per parameter would be expected to do worse with them, but that was not run.
The first fix is the reference point throughout. Singh and colleagues discuss how the start date changes NSD classification, above all when the start falls inside a migration, and a start at the range centre instead of at a random position could change the resident rates. That is a separate lesson and was left out on purpose.
The criterion is least squares on NSD, as in the curve fits it copies. The mixed migrant curve is missing, the migrant’s timing constraint is the migrateR one without its further limits on timing, and the fits come from a local search from many starting shapes, which matched or beat nls within the model on the check tracks but cannot prove that every fit is the global least squares optimum. The AR(1) repair was profiled over a grid of rho and fitted by conditional least squares, which agreed with gnls on the 2 tracks where gnls converged but is not identical to it.
Only one tau (15 days) was run for the migrants, the AR(1) refits covered three of the six values of tau, and the migrants were not refitted with AR(1) errors at all, to limit the time the page takes to build. What a floor or an AR(1) fit does to migrants with a long autocorrelation time is therefore not measured here. No true dispersers were simulated, so the disperser label is only ever an error here.
References
Bunnefeld N, Borger L, van Moorter B, Rolandsen CM, Dettki H, Solberg EJ, Ericsson G 2011 Journal of Animal Ecology 80(2):466-476 (10.1111/j.1365-2656.2010.01776.x)
Spitz DB, Hebblewhite M, Stephenson TR 2017 Ecography 40(6):788-799 (10.1111/ecog.02587)
Singh NJ, Allen AM, Ericsson G 2016 PLoS ONE 11(3):e0149594 (10.1371/journal.pone.0149594)
Cagnacci F, Focardi S, Ghisla A, van Moorter B, Merrill EH, Gurarie E, Heurich M, Mysterud A, Linnell J, Panzacchi M, May R, Nygard T, Rolandsen C, Hebblewhite M 2016 Journal of Animal Ecology 85(1):54-68 (10.1111/1365-2656.12449)
Fleming CH, Calabrese JM, Mueller T, Olson KA, Leimgruber P, Fagan WF 2014 American Naturalist 183(5):E154-E167 (10.1086/675504)