Home ranges in R: MCP versus kernel density

Estimate animal home ranges in R with minimum convex polygons and kernel density. Why the 100% MCP inflates with sample size and how bandwidth drives the KDE.

A home range is the area an animal uses over some period. From a set of relocations (GPS fixes, VHF triangulations, resightings) the question is how to turn the point cloud into an area. Two estimators dominate the older literature and still anchor most workflows: the minimum convex polygon (MCP) and the kernel density estimator (KDE). They answer slightly different questions and they fail in different ways. This post builds both from base R plus MASS, on a synthetic bilobed range, and shows the two failure modes you have to watch: the MCP grows with sample size and chases outliers, while the KDE hands the whole result to one bandwidth choice.

A synthetic range with two centres and a few excursions

Real ranges are rarely a single blob. We simulate 200 relocations from a mixture of two activity centres, then add six long excursions to stand in for the occasional foray outside the core. Coordinates are in kilometres and the data are illustrative, not a real site.

library(MASS); library(ggplot2); library(dplyr)

set.seed(4809)
n_core <- 200
comp <- rbinom(n_core, 1, 0.40)
xA <- rnorm(sum(comp == 0), 0.0, 1.20); yA <- rnorm(sum(comp == 0), 0.0, 1.20)
xB <- rnorm(sum(comp == 1), 4.0, 0.90); yB <- rnorm(sum(comp == 1), 2.5, 0.90)
core <- data.frame(x = c(xA, xB), y = c(yA, yB))

n_exc <- 6
ang <- runif(n_exc, 0, 2 * pi); rad <- runif(n_exc, 7, 10)
exc <- data.frame(x = 2 + rad * cos(ang), y = 1 + rad * sin(ang))

loc <- rbind(core, exc)
loc$kind <- c(rep("relocation", n_core), rep("excursion", n_exc))
nrow(loc)

[1] 206

Minimum convex polygon

The MCP is the smallest convex polygon that contains a chosen fraction of the relocations. The 100% version is just the convex hull of every point; chull returns the hull vertices and the shoelace formula turns them into an area.

mcp_area <- function(x, y) {
  h <- chull(x, y); xh <- x[h]; yh <- y[h]
  0.5 * abs(sum(xh * c(yh[-1], yh[1]) - c(xh[-1], xh[1]) * yh))
}
a100 <- mcp_area(loc$x, loc$y)

# 95% MCP: drop the 5% of points farthest from the centroid
cx <- mean(loc$x); cy <- mean(loc$y)
d2 <- (loc$x - cx)^2 + (loc$y - cy)^2
keep <- d2 <= quantile(d2, 0.95)
a95mcp <- mcp_area(loc$x[keep], loc$y[keep])

round(c(mcp100 = a100, mcp95 = a95mcp, ratio = a100 / a95mcp), 2)

mcp100  mcp95  ratio 
136.54  35.63   3.83 

The 100% hull covers 136.5 km squared; dropping the peripheral 5% cuts it to 35.6 km squared. A factor of almost four separates them, and the six excursions drive the gap. That sensitivity to a handful of extreme fixes is the first problem with the MCP: the estimate is defined by its most peripheral points, so one long foray can double the reported range. The 95% rule trims that, but the 5% threshold is a convention with no biological content.

The MCP grows with the number of relocations

The deeper problem is subtler. Because the convex hull can only ever expand as points are added, MCP area increases with sample size and does not settle to a stable value. We subsample the core relocations at increasing sizes and average the area over 60 random draws at each size, then do the same for a 95% kernel range for comparison.

iso_threshold <- function(z, cell, p) {
  zv <- sort(as.vector(z), decreasing = TRUE)
  zv[which(cumsum(zv) * cell >= p)[1]]
}
iso_area <- function(z, cell, p) sum(z >= iso_threshold(z, cell, p)) * cell

set.seed(7788)
ns <- seq(20, n_core, by = 20); R <- 60
xlim0 <- range(core$x) + c(-1, 1) * 2; ylim0 <- range(core$y) + c(-1, 1) * 2
acc <- lapply(ns, function(nn) {
  am <- ak <- numeric(R)
  for (r in seq_len(R)) {
    s <- core[sample(nrow(core), nn), ]
    am[r] <- mcp_area(s$x, s$y)
    kk <- kde2d(s$x, s$y, n = 80, lims = c(xlim0, ylim0))
    cc <- diff(kk$x[1:2]) * diff(kk$y[1:2])
    ak[r] <- iso_area(kk$z, cc, 0.95)
  }
  data.frame(n = nn, mcp = mean(am), kde = mean(ak))
})
acc <- do.call(rbind, acc)
round(c(mcp_n20 = acc$mcp[1], mcp_n200 = acc$mcp[10],
        kde_n20 = acc$kde[1], kde_n200 = acc$kde[10]), 1)

 mcp_n20 mcp_n200  kde_n20 kde_n200 
    19.1     38.2     53.1     45.9 

# Second arm, on exactly the same subsamples: re-seeding replays the draws
# above, but now the bandwidth is held at the full-sample value instead of
# being recomputed inside each subsample.
h_full <- c(bandwidth.nrd(core$x), bandwidth.nrd(core$y))
set.seed(7788)
fixed_bw <- lapply(ns, function(nn) {
  af <- hx <- numeric(R)
  for (r in seq_len(R)) {
    s <- core[sample(nrow(core), nn), ]
    hx[r] <- bandwidth.nrd(s$x)
    kk <- kde2d(s$x, s$y, n = 80, lims = c(xlim0, ylim0), h = h_full)
    cc <- diff(kk$x[1:2]) * diff(kk$y[1:2])
    af[r] <- iso_area(kk$z, cc, 0.95)
  }
  data.frame(n = nn, kde_fix = mean(af), h_x = mean(hx))
})
fixed_bw <- do.call(rbind, fixed_bw)
acc$kde_fix <- fixed_bw$kde_fix; acc$h_x <- fixed_bw$h_x

data.frame(n = acc$n, kde = round(acc$kde, 1),
           kde_fix = round(acc$kde_fix, 1), h_x = round(acc$h_x, 2))

     n  kde kde_fix  h_x
1   20 53.1    37.9 4.71
2   40 53.1    41.6 4.36
3   60 51.7    43.8 4.02
4   80 50.5    44.9 3.78
5  100 48.8    44.8 3.62
6  120 48.2    45.3 3.50
7  140 47.5    45.6 3.37
8  160 46.8    45.6 3.29
9  180 46.3    45.8 3.22
10 200 45.9    45.9 3.15

round(c(h_ratio = acc$h_x[10] / acc$h_x[1], nrd_rule = 10^(-1 / 5)), 3)

 h_ratio nrd_rule 
   0.669    0.631 

From 20 to 200 relocations the mean MCP area climbs from 19.1 to 38.2 km squared, a doubling with no sign of levelling off. The kernel range moves the other way, from 53.1 to 45.9 km squared, but that fall belongs to the bandwidth rule rather than to the estimator. bandwidth.nrd scales as the sample size to the power minus one fifth, so it shrinks as fixes accumulate: on the x axis it drops from 4.71 at 20 relocations to 3.15 at 200, a ratio of 0.669 against the 0.631 that rule predicts for a ten-fold change in n. Hold the bandwidth at the full-sample value and the kernel area rises with sample size as well, from 37.9 to 45.9 km squared, flattening out beyond about 80 relocations while the hull keeps climbing.

The contrast that survives is narrower than it first looks: the kernel estimate settles where the MCP does not, but it settles because the selector trades bandwidth against sample size, not because the estimator is indifferent to effort. This is why MCP areas from studies with different tracking effort are not comparable: more fixes mean a larger reported range, all else equal. Reviews of home-range methods have made this point repeatedly (Harris 1990; Powell 2000).

lev <- c("100% MCP", "95% KDE, bandwidth per subsample", "95% KDE, fixed bandwidth")
acc_long <- rbind(
  data.frame(n = acc$n, area = acc$mcp, estimator = lev[1]),
  data.frame(n = acc$n, area = acc$kde, estimator = lev[2]),
  data.frame(n = acc$n, area = acc$kde_fix, estimator = lev[3]))
acc_long$estimator <- factor(acc_long$estimator, levels = lev)
ggplot(acc_long, aes(n, area, colour = estimator)) +
  geom_line(linewidth = 0.8) + geom_point(size = 1.8) +
  scale_colour_manual(
    values = setNames(c(te_red, te_forest, te_gold), lev), name = NULL) +
  labs(title = "Home-range area against number of relocations",
       x = "number of relocations", y = expression(area ~ (km^2))) +
  theme_te(13)

Figure 1: Mean home-range area against the number of relocations, averaged over 60 random subsamples. The 100% MCP keeps climbing; the kernel range falls only while the bandwidth is recomputed inside each subsample, and rises when it is held at the full-sample value.

Kernel density and volume-based isopleths

The KDE replaces the hard hull with a smooth utilisation distribution: a density surface over the plane, estimated by placing a kernel on each relocation and summing. MASS::kde2d does this on a grid, and by default it picks a bandwidth from the normal reference rule for each axis.

xr <- range(loc$x) + c(-1, 1) * 2; yr <- range(loc$y) + c(-1, 1) * 2
href <- c(bandwidth.nrd(loc$x), bandwidth.nrd(loc$y))
k <- kde2d(loc$x, loc$y, n = 120, lims = c(xr, yr), h = href)
cell <- diff(k$x[1:2]) * diff(k$y[1:2])
round(c(h_x = href[1], h_y = href[2], integral = sum(k$z) * cell), 3)

     h_x      h_y integral 
   3.429    2.614    1.000 

To get a home range from a density surface you take an isopleth: the contour enclosing a given share of the total volume. The 95% isopleth is the smallest region holding 95% of the utilisation, and the 50% isopleth marks the core. Sort the grid densities from high to low, accumulate their volume, and read off the density level where the running total crosses the target.

thr95 <- iso_threshold(k$z, cell, 0.95); thr50 <- iso_threshold(k$z, cell, 0.50)
a95kde <- iso_area(k$z, cell, 0.95); a50kde <- iso_area(k$z, cell, 0.50)
round(c(kde95 = a95kde, kde50 = a50kde), 2)

kde95 kde50 
59.50 13.86 

The 95% kernel range is 59.5 km squared and the 50% core is 13.9 km squared. Unlike the convex hull, the kernel surface can represent the gap between the two centres and the separate cores, which a single polygon cannot.

hull_idx <- chull(loc$x, loc$y)
hull_df  <- loc[c(hull_idx, hull_idx[1]), c("x", "y")]
grid_df  <- expand.grid(x = k$x, y = k$y); grid_df$z <- as.vector(k$z)

ggplot() +
  geom_point(data = loc, aes(x, y, shape = kind, colour = kind), size = 1.6, alpha = 0.8) +
  geom_path(data = hull_df, aes(x, y), colour = te_red, linewidth = 0.7) +
  geom_contour(data = grid_df, aes(x, y, z = z), breaks = thr95, colour = te_forest, linewidth = 0.8) +
  geom_contour(data = grid_df, aes(x, y, z = z), breaks = thr50, colour = te_gold, linewidth = 0.8) +
  scale_colour_manual(values = c(relocation = te_faint, excursion = te_red), name = NULL) +
  scale_shape_manual(values = c(relocation = 16, excursion = 4), name = NULL) +
  coord_equal() +
  labs(title = "Home range: 100% MCP versus kernel isopleths",
       x = "easting (km)", y = "northing (km)") +
  theme_te(13)

Figure 2: Relocations with the 100% MCP hull (red), the 95% kernel isopleth (green) and the 50% core (gold). Crosses mark the six excursions that stretch the hull.

Bandwidth is the choice that matters

The kernel result is only as good as the bandwidth. The normal reference rule assumes a single roughly Gaussian blob, so for a multimodal range it tends to oversmooth: the surface spreads across the gap and the isopleth inflates. Scaling the reference bandwidth up and down shows how much rides on it.

bw_area <- function(mult) {
  kk <- kde2d(loc$x, loc$y, n = 120, lims = c(xr, yr), h = href * mult)
  cc <- diff(kk$x[1:2]) * diff(kk$y[1:2])
  iso_area(kk$z, cc, 0.95)
}
round(sapply(c(0.6, 1.0, 1.5), bw_area), 2)

[1] 45.15 59.50 81.28

At 0.6 times the reference the 95% range is 45.2 km squared; at the reference it is 59.5; at 1.5 times it is 81.3. The estimate almost doubles across a plausible band of smoothing, with the data unchanged. This is the central trade-off: too small a bandwidth breaks the range into islands around individual fixes, too large a bandwidth smears it into one oversized blob. Least-squares cross-validation and plug-in selectors try to choose objectively, and their behaviour has been studied at length (Worton 1989; Seaman and Powell 1996), but no rule removes the judgement entirely.

Duration sets the precision, not the fix count

Everything so far holds the relocations fixed and asks what an estimator makes of them. The schedule that produced them has a failure mode of its own, and it is the expensive one: collar battery and memory go on the fix interval, while the answer mostly depends on how long the collar stays on the animal. Consecutive positions are correlated, and the time constant of that correlation, tau, is roughly how long the animal takes to forget where it was. A fix taken well inside tau repeats most of what the previous one already said.

A two-dimensional Ornstein-Uhlenbeck process makes the argument checkable. The animal is pulled back towards a centre with time constant tau, and its stationary distribution is a bivariate normal with standard deviation sigma on each axis. Positions separated by an interval dt satisfy an exact update, x_next = mu + exp(-dt / tau) * (x - mu) + sigma * sqrt(1 - exp(-2 * dt / tau)) * z, so the simulated track carries no discretisation error at any fix interval, and one rule generates the schedule whether the collar fires every minute or once a week.

q95 <- qchisq(0.95, 2)
tau_h <- 6          # position autocorrelation time, hours
sigma_km <- 1.2     # stationary standard deviation per axis, km
area_true <- pi * q95 * sigma_km^2
range_diam <- 2 * sqrt(q95) * sigma_km

# exact OU update, run as a recursive filter on the deviations from the centre
ou_dev <- function(n, dt, tau, sigma) {
  aa <- exp(-dt / tau)
  as.numeric(stats::filter(rnorm(n, 0, sigma * sqrt(1 - aa^2)), aa,
                           method = "recursive", init = rnorm(1, 0, sigma)))
}
# the model's own target, estimated its own way: isotropic normal 95% area
area_norm <- function(x, y) {
  vx <- mean((x - mean(x))^2); vy <- mean((y - mean(y))^2)
  pi * q95 * (vx + vy) / 2
}
# expected area over true area for n fixes dt apart, from the discrete
# autocovariance rather than its long-track limit
ratio_exact <- function(n, dt, tau) {
  rho <- exp(-dt / tau); kk <- seq_len(n - 1)
  1 - (n + 2 * sum((n - kk) * rho^kk)) / n^2
}

set.seed(60481)
chk <- ou_dev(2e5, 0.25 * tau_h, tau_h, sigma_km)
chk_sd <- sd(chk); chk_acf <- cor(chk[-1], chk[-length(chk)])
round(c(sd_sim = chk_sd, sd_target = sigma_km, acf1_sim = chk_acf,
        acf1_target = exp(-0.25), area_true = area_true), 4)

     sd_sim   sd_target    acf1_sim acf1_target   area_true 
     1.1917      1.2000      0.7770      0.7788     27.1047 

The simulator returns the spread and the lag one correlation it was asked for, so the machinery is right. With tau at 6 hours and sigma at 1.2 km, the true 95% area is 27.10 km squared and the range is 5.87 km across.

Because the movement model is known, the target is known too, and the estimator can be matched to it: average the two coordinate variances taken about the sample mean, which is the isotropic form of the normal 95% ellipse. The full covariance version adds a cross term carrying a bias of its own, and this model is isotropic, so the simple form is the one to test theory against. What comes out is a statement about the schedule and not about the estimator: a short track has seen a small part of the range whatever you fit to it afterwards.

Two consequences follow in closed form. The first is bias. The mean of an OU track has variance 2 * sigma^2 * tau / T once the track is long against tau, and the variance taken about that estimated mean falls short by exactly that much, so the expected area is 1 - 2 * tau / T of the truth. The second is precision. Squared deviations of a Gaussian process decorrelate twice as fast as the positions do, which puts the coefficient of variation of one coordinate variance at sqrt(2 * tau / T). The area averages the two coordinates, and under an isotropic model those two variances are independent, so averaging them halves the variance once more and leaves the coefficient of variation of the area at sqrt(tau / T). Neither expression mentions the number of fixes.

The first sweep holds the interval at a quarter of tau and lengthens the track.

boot_cv <- function(v, B = 400) {
  bs <- replicate(B, {u <- sample(v, length(v), TRUE); sd(u) / mean(u)})
  sd(bs)
}

set.seed(31607)
dt_base <- 0.25; R_rep <- 1000
tmult <- c(2, 5, 10, 25, 50, 100, 400)
dur <- lapply(tmult, function(tm) {
  nfix <- tm / dt_base
  a <- replicate(R_rep, area_norm(ou_dev(nfix, dt_base * tau_h, tau_h, sigma_km),
                                  ou_dev(nfix, dt_base * tau_h, tau_h, sigma_km)))
  rr <- a / area_true; cv <- sd(rr) / mean(rr)
  data.frame(t_tau = tm, n_fix = nfix, ratio = mean(rr), se = sd(rr) / sqrt(R_rep),
             asym = 1 - 2 / tm, exact = ratio_exact(nfix, dt_base * tau_h, tau_h),
             cv = cv, cv_se = boot_cv(rr), cv_asym = sqrt(1 / tm))
})
dur <- do.call(rbind, dur)
cv_gain <- dur$cv[6] / dur$cv[5]
cv_gain_se <- cv_gain * sqrt((dur$cv_se[6] / dur$cv[6])^2 + (dur$cv_se[5] / dur$cv[5])^2)
round(dur, 4)

  t_tau n_fix  ratio     se  asym  exact     cv  cv_se cv_asym
1     2     8 0.4244 0.0075 0.000 0.4249 0.5626 0.0287  0.7071
2     5    20 0.6814 0.0092 0.600 0.6770 0.4275 0.0152  0.4472
3    10    40 0.8160 0.0078 0.800 0.8189 0.3041 0.0089  0.3162
4    25   100 0.9263 0.0062 0.920 0.9228 0.2133 0.0042  0.2000
5    50   200 0.9617 0.0044 0.960 0.9606 0.1447 0.0033  0.1414
6   100   400 0.9837 0.0032 0.980 0.9801 0.1040 0.0023  0.1000
7   400  1600 0.9948 0.0016 0.995 0.9950 0.0508 0.0011  0.0500

The design runs from 2 tau to 400 tau of tracking, 1000 replicate tracks each, and both closed forms hold over most of it. At 25 tau the mean area is 0.926 of the truth with a Monte Carlo standard error of 0.006, against 0.920 from 1 - 2 * tau / T; the coefficient of variation is 0.213 against 0.200 from sqrt(tau / T). At 400 tau the same pair reads 0.995 against 0.995 and 0.0508 against 0.0500. The asymptotic bias formula breaks down at the short end, where it has to: at 2 tau it predicts an area of 0.0, the simulation returns 0.424, and the exact discrete expression returns 0.425. The square root law is the part worth carrying around: doubling the duration from 50 tau to 100 tau multiplies the coefficient of variation by 0.719, standard error 0.023, against the 0.707 that halving the variance would give.

The second sweep holds the duration at 25 tau and moves the fix interval instead, from coarser than tau down to a sixteenth of it. The 100% MCP from the earlier section rides along on the same tracks.

set.seed(90412)
dtm <- c(2.5, 1.25, 0.5, 0.25, 0.125, 0.0625); t_hold <- 25
pace <- lapply(dtm, function(dm) {
  nfix <- t_hold / dm; an <- am <- numeric(R_rep)
  for (r in seq_len(R_rep)) {
    xx <- ou_dev(nfix, dm * tau_h, tau_h, sigma_km)
    yy <- ou_dev(nfix, dm * tau_h, tau_h, sigma_km)
    an[r] <- area_norm(xx, yy); am[r] <- mcp_area(xx, yy)
  }
  rr <- an / area_true; cv <- sd(rr) / mean(rr)
  data.frame(dt_tau = dm, n_fix = nfix, ratio = mean(rr), se = sd(rr) / sqrt(R_rep),
             exact = ratio_exact(nfix, dm * tau_h, tau_h), cv = cv, cv_se = boot_cv(rr),
             mcp = mean(am), mcp_se = sd(am) / sqrt(R_rep))
})
pace <- do.call(rbind, pace)
fix_gain <- pace$cv[6] / pace$cv[5]
fix_gain_se <- fix_gain * sqrt((pace$cv_se[6] / pace$cv[6])^2 + (pace$cv_se[5] / pace$cv[5])^2)
mcp_gain <- pace$mcp[6] / pace$mcp[4]
mcp_gain_se <- mcp_gain * sqrt((pace$mcp_se[6] / pace$mcp[6])^2 + (pace$mcp_se[4] / pace$mcp[4])^2)
n_eff <- 1 / pace$cv[6]^2
z_share <- (pace$ratio[4] - dur$ratio[4]) / sqrt(pace$se[4]^2 + dur$se[4]^2)

# what T / tau independent relocations would give, for comparison
set.seed(41137)
iid <- replicate(R_rep, area_norm(rnorm(t_hold, 0, sigma_km), rnorm(t_hold, 0, sigma_km)))
iid_cv <- sd(iid) / mean(iid)
round(pace, 4)

  dt_tau n_fix  ratio     se  exact     cv  cv_se     mcp mcp_se
1 2.5000    10 0.8853 0.0098 0.8841 0.3509 0.0101  7.5541 0.0937
2 1.2500    20 0.9261 0.0072 0.9127 0.2475 0.0062 12.0620 0.1023
3 0.5000    50 0.9154 0.0059 0.9215 0.2047 0.0050 16.9432 0.1128
4 0.2500   100 0.9242 0.0060 0.9228 0.2069 0.0047 20.0774 0.1216
5 0.1250   200 0.9169 0.0058 0.9231 0.2004 0.0055 22.4738 0.1346
6 0.0625   400 0.9285 0.0060 0.9232 0.2032 0.0045 24.6120 0.1347

The first step down is worth paying for, because sampling more coarsely than tau does throw information away: going from an interval of 2.50 tau to 1.25 tau drops the coefficient of variation from 0.351 to 0.247. After that the dial stops working. Over the last eightfold increase in fixes, from 50 to 400 on the same length of track, the coefficient of variation reads 0.205, 0.207, 0.200 and 0.203, and the final doubling of the fix rate multiplies it by 1.014 with a standard error of 0.036, which is one. The mean area is flat across the same stretch, and the exact prediction is flat with it, at 0.921 and 0.923. The two sweeps share the design point at 25 tau with an interval of 0.25 tau, measured twice under independent seeds, and the two values sit 0.24 standard errors apart.

That flat column has a plain reading. A coefficient of variation of 0.203 off 400 fixes is what 24.2 independent relocations would have delivered, and 25 independent relocations do deliver 0.202 on the same estimator. The number of independent looks at the range is the duration in units of tau, 25 here, and it takes no notice of how many rows the collar wrote. The MCP behaves as the earlier section said it would: at a fixed duration, quadrupling the fix rate from 100 to 400 fixes raises the mean hull area by a factor of 1.226, standard error 0.010, on tracking that contains no more information about the range than it did before.

arm <- c("more fixes from a longer track", "more fixes from a faster fix rate")
qn <- c("area estimate / true area", "coefficient of variation")
d2 <- dur[dur$t_tau >= 5, ]
sched <- rbind(
  data.frame(n_fix = d2$n_fix, value = d2$ratio, arm = arm[1], quantity = qn[1], src = "simulation"),
  data.frame(n_fix = d2$n_fix, value = d2$cv, arm = arm[1], quantity = qn[2], src = "simulation"),
  data.frame(n_fix = pace$n_fix, value = pace$ratio, arm = arm[2], quantity = qn[1], src = "simulation"),
  data.frame(n_fix = pace$n_fix, value = pace$cv, arm = arm[2], quantity = qn[2], src = "simulation"),
  data.frame(n_fix = d2$n_fix, value = d2$asym, arm = arm[1], quantity = qn[1], src = "closed form"),
  data.frame(n_fix = d2$n_fix, value = d2$cv_asym, arm = arm[1], quantity = qn[2], src = "closed form"),
  data.frame(n_fix = pace$n_fix, value = 1 - 2 / t_hold, arm = arm[2], quantity = qn[1], src = "closed form"),
  data.frame(n_fix = pace$n_fix, value = sqrt(1 / t_hold), arm = arm[2], quantity = qn[2], src = "closed form"))
sched$arm <- factor(sched$arm, levels = arm)
sched$quantity <- factor(sched$quantity, levels = qn)

ggplot(sched, aes(n_fix, value, colour = arm, linetype = src)) +
  geom_line(linewidth = 0.8) +
  geom_point(data = subset(sched, src == "simulation"), size = 1.8, show.legend = FALSE) +
  facet_wrap(~ quantity, scales = "free_y") +
  scale_x_log10() +
  scale_colour_manual(values = setNames(c(te_forest, te_gold), arm), name = NULL) +
  scale_linetype_manual(values = c(simulation = "solid", `closed form` = "22"), name = NULL) +
  guides(colour = guide_legend(order = 1), linetype = guide_legend(order = 2)) +
  labs(title = "What extra fixes buy, and what they do not",
       x = "number of fixes (log scale)", y = "dimensionless ratio") +
  theme_te(13) +
  theme(legend.position = "bottom", legend.box = "vertical",
        legend.margin = margin(0, 0, 0, 0))

Figure 3: Two ways of buying more fixes, plotted against the number of fixes on a log scale. Green: a longer track at a fixed interval of a quarter of tau. Gold: a faster fix rate over a fixed 25 tau of tracking. Dashed lines are the closed forms. The shortest track in the table is left out because the asymptotic bias formula is meaningless there.

This sits alongside the accumulation section rather than against it. That section resampled independent relocations and let their number grow, so its x axis was already an independent-fix count; this one holds the window fixed and lets the fixes crowd together inside it. Both find the hull inflating whenever the fix count rises, and they differ only in what a fix is worth: an independent relocation adds information, a fix taken inside tau mostly does not, so on a real collar the quantity that belongs on that earlier x axis is the duration in units of tau, not the number of rows in the file.

Now the hedge, and it is not a small one. All of the above assumes the range sits still. Let the centre drift and the recommendation reverses, because a longer track then answers a different question: the ground covered over a season rather than the area in use at any point within it. The same OU model gets a centre sliding along the x axis as an amplitude times sin(2 * pi * t / P), the target stays the within-season area, and the score is relative root mean squared error against it.

set.seed(55231)
amp_km <- 1.5 * sigma_km; per_h <- 250 * tau_h; R_sea <- 400
tm_sea <- c(5, 25, 100, 500, 2000)
sea <- lapply(tm_sea, function(tm) {
  nfix <- tm / dt_base; dth <- dt_base * tau_h; tt <- seq_len(nfix) * dth
  ef <- es <- numeric(R_sea)
  for (r in seq_len(R_sea)) {
    xx <- ou_dev(nfix, dth, tau_h, sigma_km); yy <- ou_dev(nfix, dth, tau_h, sigma_km)
    ef[r] <- area_norm(xx, yy) / area_true - 1
    es[r] <- area_norm(xx + amp_km * sin(2 * pi * tt / per_h + runif(1, 0, 2 * pi)),
                       yy) / area_true - 1
  }
  data.frame(t_tau = tm, rmse_fix = sqrt(mean(ef^2)),
             se_fix = sd(ef^2) / (2 * sqrt(mean(ef^2)) * sqrt(R_sea)),
             rmse_sea = sqrt(mean(es^2)),
             se_sea = sd(es^2) / (2 * sqrt(mean(es^2)) * sqrt(R_sea)))
})
sea <- do.call(rbind, sea)
infl_pred <- amp_km^2 / (4 * sigma_km^2)
sea_worse <- sea$rmse_sea[5] / sea$rmse_fix[5]
swing_frac <- 2 * amp_km / range_diam
per_days <- per_h / 24
round(sea, 4)

  t_tau rmse_fix se_fix rmse_sea se_sea
1     5   0.4242 0.0082   0.4225 0.0081
2    25   0.2000 0.0058   0.1994 0.0059
3   100   0.1033 0.0033   0.3052 0.0095
4   500   0.0464 0.0016   0.5620 0.0041
5  2000   0.0238 0.0008   0.5601 0.0019

The shift has an amplitude of 1.80 km and a period of 62.5 days, which is 250 tau; the centre swings back and forth over 61 per cent of the range diameter, so this is a seasonal drift and not a migration. Each of the 400 replicate tracks per design starts at a random phase of the season. With the centre held still the error falls all the way along the grid, from 0.424 at 5 tau to 0.0238 at 2000 tau. With the centre drifting the two are indistinguishable up to 25 tau, 0.199 against 0.200 with standard errors of 0.006 and 0.006, and after that the drift takes over: 0.305 at 100 tau, 0.562 at 500 tau and 0.560 at 2000 tau, which is 23 times worse than a fixed-centre track of the same length. That last value is bias and not noise: a centre sweeping a sine wave adds a quarter of its squared amplitude to the average of the two coordinate variances, which predicts an inflation of 0.5625 against the measured 0.5601, standard error 0.0019.

lab <- c("fixed range centre", "seasonal shift of the centre")
sea_long <- rbind(
  data.frame(t_tau = sea$t_tau, rmse = sea$rmse_fix, se = sea$se_fix, model = lab[1]),
  data.frame(t_tau = sea$t_tau, rmse = sea$rmse_sea, se = sea$se_sea, model = lab[2]))
sea_long$model <- factor(sea_long$model, levels = lab)

ggplot(sea_long, aes(t_tau, rmse, colour = model)) +
  geom_vline(xintercept = per_h / tau_h, colour = te_faint, linetype = "22", linewidth = 0.4) +
  geom_line(linewidth = 0.8) + geom_point(size = 1.9) +
  geom_errorbar(aes(ymin = rmse - se, ymax = rmse + se), width = 0.06, linewidth = 0.5) +
  scale_x_log10() +
  scale_colour_manual(values = setNames(c(te_forest, te_red), lab), name = NULL) +
  labs(title = "Longer tracks stop paying when the centre moves",
       x = expression("tracking duration in units of" ~ tau ~ "(log scale)"),
       y = "relative RMSE of the area estimate") +
  theme_te(13) +
  theme(legend.position = "bottom")

Figure 4: Relative RMSE of the area estimate against tracking duration, with a fixed range centre and with a centre that slides seasonally. The vertical line marks one seasonal period. The two curves lie on top of each other until the track is long enough to feel the drift.

So the advice is not to track for longer. It is to track for longer than the autocorrelation time and for less than whatever moves the range, and to say which period the estimate refers to. Of the durations on this grid the best under seasonal drift is 25 tau, long enough for the sampling error to have come down and short enough for the drift not to have arrived; where that optimum falls depends on the amplitude and the period, and a wider swing moves it to the left. If the range shifts across the year, the honest product is one estimate per season with the duration reported alongside it, not one estimate per collar.

Which to use

The MCP is quick, needs no tuning and is still the standard for a crude outer boundary, but its area depends on sample size and on the most extreme fixes, so it is a poor choice for comparing ranges across animals or studies with uneven effort. The KDE gives a proper utilisation distribution with a defensible core, handles multimodal ranges, and is far less sensitive to sample size, partly because the reference bandwidth rule shrinks as fixes accumulate and offsets the rise that remains at a fixed bandwidth, at the cost of a bandwidth decision that changes the answer. Report which estimator and, for the kernel, which bandwidth selector you used; without that, a home-range area is not reproducible. Both estimators also assume the relocations are an unbiased sample of use, which autocorrelated tracking data and gappy fix schedules can violate.

References

Mohr 1947 American Midland Naturalist 37(1):223-249 (10.2307/2421652)

Worton 1989 Ecology 70(1):164-168 (10.2307/1938423)

Harris, Cresswell, Forde, Trewhella, Woollard & Wray 1990 Mammal Review 20(2-3):97-123 (10.1111/j.1365-2907.1990.tb00106.x)

Seaman & Powell 1996 Ecology 77(7):2075-2085 (10.2307/2265701)

Powell 2000, in Research Techniques in Animal Ecology (Boitani & Fuller, eds), Columbia University Press, ISBN 978-0-231-11341-2

Silverman 1986 Density Estimation for Statistics and Data Analysis, Chapman and Hall, ISBN 978-0-412-24620-3

Related tutorials

• Complete spatial randomness and quadrat tests
• Step lengths and turning angles
• Correlated random walks and net displacement
• Isotopic niche width: hulls and ellipses