Staying time and the REST model in R

camera traps

abundance

movement ecology

ecology tutorial

The REST model estimates animal density from how long animals stay in view, with no external speed. Little’s law and Cauchy’s chord identity, tested in R.

Author

Tidy Ecology

Published

2026-07-16

The random encounter model has one weak joint: it needs the average speed, and the average has to include every hour the animal spends lying still. You cannot get that from the camera without an assumption about activity, and you cannot get it from a GPS collar without collaring something.

The random encounter and staying time model (Nakashima et al. 2018) removes that joint. It uses the same cameras and the same triggers, but instead of asking how fast animals cross the detection zone it asks how long they stay in it. The estimator is:

\[D = \frac{y}{t} \cdot \frac{\mathbb{E}[T_{\text{stay}}]}{a}\]

where $a$ is the area of the detection zone and $T_{\text{stay}}$ is the time an animal spends in view per visit. No speed. No activity level. Both are measured for you, implicitly, by the stopwatch.

Why that works: Little’s law

The formula is not a camera trap result at all. It is Little’s law, which says that for any stable system, the average number of customers inside equals the arrival rate times the average time each one spends there. Queueing theorists have used it since 1961 to reason about supermarket tills.

Point it at a detection zone. The average number of animals inside is $D a$, because the gas is uniform. The arrival rate is $y/t$. The average time each one spends there is $\mathbb{E}[T_{\text{stay}}]$. So

\[D a = \frac{y}{t}\,\mathbb{E}[T_{\text{stay}}]\]

and that is the whole derivation. Notice what is not in it: no assumption about how the animals move, no speed, no direction distribution, not even the shape of the zone. Little’s law needs stationarity and nothing else.

Set it beside the random encounter model, which says $y/t = D v r(2+\theta)/\pi$, and one substitution gives a prediction about staying time:

\[\mathbb{E}[T_{\text{stay}}] = \frac{\pi a}{v L}, \qquad L = r(2+\theta)\]

That is Cauchy’s mean chord formula wearing a wildlife hat: the average chord an isotropic straight line cuts across a convex region has length $\pi a / L$, and dividing by speed turns a length into a time. So the two models are the same equation read in two directions. The random encounter model measures $v$ elsewhere and predicts the count; the staying time model measures the count and the clock, and lets them tell it about $v$.

The rig

Same gas as the encounter rate post: independent walkers on a torus, exact boundary crossings, no time-step discretisation to apologise for. The extra piece here is the stopwatch.

inside <- function(P, r, half) {
  (P[, 1]^2 + P[, 2]^2 <= r^2) & (abs(atan2(P[, 2], P[, 1])) <= half)
}

crossings <- function(P1, D, r, half) {
  m <- nrow(P1)
  U <- matrix(Inf, m, 4); IN <- matrix(FALSE, m, 4)
  a <- D[, 1]^2 + D[, 2]^2
  b <- 2 * (P1[, 1] * D[, 1] + P1[, 2] * D[, 2])
  cc <- P1[, 1]^2 + P1[, 2]^2 - r^2
  disc <- b^2 - 4 * a * cc
  ok <- disc > 0 & a > 0
  sq <- ifelse(ok, sqrt(pmax(disc, 0)), NA_real_)
  for (k in 1:2) {                                        # the arc
    u <- (-b + c(-1, 1)[k] * sq) / (2 * a)
    Q <- cbind(P1[, 1] + u * D[, 1], P1[, 2] + u * D[, 2])
    good <- ok & !is.na(u) & u > 0 & u <= 1 & abs(atan2(Q[, 2], Q[, 1])) <= half
    good[is.na(good)] <- FALSE
    U[good, k] <- u[good]
    IN[good, k] <- ((D[, 1] * Q[, 1] + D[, 2] * Q[, 2]) < 0)[good]
  }
  for (k in 1:2) {                                        # the two radii
    s <- c(1, -1)[k] * half
    e <- c(cos(s), sin(s))
    nrm <- c(cos(s + c(1, -1)[k] * pi / 2), sin(s + c(1, -1)[k] * pi / 2))
    cD <- D[, 1] * e[2] - D[, 2] * e[1]
    u <- -(P1[, 1] * e[2] - P1[, 2] * e[1]) / cD
    tt <- (P1[, 1] + u * D[, 1]) * e[1] + (P1[, 2] + u * D[, 2]) * e[2]
    good <- abs(cD) > 1e-12 & !is.na(u) & u > 0 & u <= 1 & tt >= 0 & tt <= r
    good[is.na(good)] <- FALSE
    U[good, 2 + k] <- u[good]
    IN[good, 2 + k] <- ((D[, 1] * nrm[1] + D[, 2] * nrm[2]) < 0)[good]
  }
  list(U = U, IN = IN)
}

An animal is active or resting, switching between the two as a Markov chain on the step grid. Resting means the step length is zero, which the crossing routine handles for free: a zero-length segment cannot cross anything. The proportion of time active follows from the two switching probabilities.

walk_one <- function(nstep, dt, v, W, r, half, rho, p_ar = 0, p_ra = 1) {
  ## active / resting state sequence, built bout by bout (p_ar = 0 means never rests)
  st <- if (p_ar <= 0) rep(1L, nstep) else {
    z <- integer(0); s <- 1L
    while (length(z) < nstep) {
      z <- c(z, rep(s, rgeom(1, if (s == 1L) p_ar else p_ra) + 1L)); s <- 1L - s
    }
    z[seq_len(nstep)]
  }

  u <- runif(nstep)
  turn <- 2 * atan(((1 - rho) / (1 + rho)) * tan(pi * (u - 0.5)))
  phi <- runif(1, 0, 2 * pi) + cumsum(turn)
  step <- v * dt * st                                     # zero while resting
  D <- cbind(step * cos(phi), step * sin(phi))
  p0 <- c(runif(1, -W / 2, W / 2), runif(1, -W / 2, W / 2))
  X <- p0[1] + cumsum(D[, 1]); Y <- p0[2] + cumsum(D[, 2])
  P1 <- cbind(c(p0[1], X[-nstep]), c(p0[2], Y[-nstep]))
  R1 <- cbind(((P1[, 1] + W / 2) %% W) - W / 2, ((P1[, 2] + W / 2) %% W) - W / 2)

  frac_in <- mean(inside(R1, r, half))                    # time-average occupancy
  idx <- which((R1[, 1]^2 + R1[, 2]^2) < (r + v * dt)^2)
  if (!length(idx)) return(list(visits = numeric(0), frac_in = frac_in, active = mean(st)))
  cr <- crossings(R1[idx, , drop = FALSE], D[idx, , drop = FALSE], r, half)
  hit <- which(is.finite(cr$U), arr.ind = TRUE)
  if (!nrow(hit)) return(list(visits = numeric(0), frac_in = frac_in, active = mean(st)))
  tt <- (idx[hit[, 1]] - 1 + cr$U[hit]) * dt
  o <- order(tt); tt <- tt[o]; inw <- cr$IN[hit][o]
  t_in <- tt[inw]; t_out <- tt[!inw]
  if (!length(t_in)) return(list(visits = numeric(0), frac_in = frac_in, active = mean(st)))
  j <- findInterval(t_in, t_out) + 1L                     # the exit that closes each entry
  ok <- j <= length(t_out)                                # drop a visit still open at the end
  list(visits = t_out[j[ok]] - t_in[ok], frac_in = frac_in, active = mean(st))
}

gas <- function(n, nstep, dt, v, W, r, half, rho, p_ar = 0, p_ra = 1) {
  out <- lapply(seq_len(n), function(i) walk_one(nstep, dt, v, W, r, half, rho, p_ar, p_ra))
  list(visits = unlist(lapply(out, `[[`, "visits")),
       dur = vapply(out, function(z) sum(z$visits), numeric(1)),    # per animal
       nvis = vapply(out, function(z) length(z$visits), numeric(1)),
       frac_in = vapply(out, `[[`, numeric(1), "frac_in"),
       active = vapply(out, `[[`, numeric(1), "active"))
}

## the mean staying time is a ratio of two per-animal totals, so its standard error
## comes from resampling animals, which also allows for any clustering of visits
boot_stay <- function(g, B = 600) {
  n <- length(g$dur)
  vapply(seq_len(B), function(b) {
    i <- sample.int(n, n, replace = TRUE)
    sum(g$dur[i]) / sum(g$nvis[i])
  }, numeric(1))
}

Four identities, one simulation

Start with a gas that never rests. There are four claims stacked on top of each other and the simulation checks each one separately: that the gas really is uniform, that Little’s law holds, that the mean staying time matches Cauchy’s chord formula, and that the resulting density estimate is right.

W <- 20; r <- 4; half <- 20 * pi / 180; v <- 1; dt <- 0.5
A <- W^2; a_zone <- half * r^2; L <- 2 * r + 2 * r * half
nanim <- 120; nstep <- 20000; Tt <- nstep * dt; Dtrue <- nanim / A

set.seed(2701)
g <- gas(nanim, nstep, dt, v, W, r, half, rho = 0.9)
lambda <- length(g$visits) / Tt                    # arrivals per second, whole population
Wbar <- mean(g$visits)                             # mean staying time
Lnum <- sum(g$frac_in)                             # mean number inside, from positions

data.frame(
  claim = c("gas is uniform: mean number inside = D a",
            "Little's law: arrival rate x staying time = mean number inside",
            "Cauchy chord: mean staying time = pi a / (v L)",
            "REST estimate of density"),
  measured = c(Lnum, lambda * Wbar, Wbar, lambda * Wbar / a_zone),
  predicted = c(Dtrue * a_zone, Lnum, pi * a_zone / (v * L), Dtrue))

                                                           claim  measured
1                       gas is uniform: mean number inside = D a 1.6662000
2 Little's law: arrival rate x staying time = mean number inside 1.6673568
3                 Cauchy chord: mean staying time = pi a / (v L) 1.6258964
4                                       REST estimate of density 0.2985391
  predicted
1  1.675516
2  1.666200
3  1.625751
4  0.300000

Every row lands. The third is the one worth pausing on: the mean staying time is 1.626 seconds against a predicted 1.626, and that prediction came from a nineteenth-century integral geometry result about random lines crossing a convex set.

Now vary the tortuosity. The encounter rate post showed that path persistence does not change how often animals arrive. Little’s law then forces something less obvious: it cannot change how long they stay either, because the occupancy and the arrival rate are both fixed, and the staying time is their ratio.

set.seed(2702)
tor <- do.call(rbind, lapply(c(0, 0.5, 0.9, 0.98), function(rh) {
  gg <- gas(150, nstep, dt, v, W, r, half, rho = rh)
  bs <- boot_stay(gg)
  data.frame(rho = rh, n_visits = length(gg$visits), mean_stay = mean(gg$visits),
             boot_se = sd(bs), CV_stay = sd(gg$visits) / mean(gg$visits))
}))
tor$predicted <- pi * a_zone / (v * L)
tor$ratio <- tor$mean_stay / tor$predicted
tor$z <- (tor$mean_stay - tor$predicted) / tor$boot_se
tor

   rho n_visits mean_stay     boot_se   CV_stay predicted     ratio           z
1 0.00    13057  1.605594 0.019604282 1.2827629  1.625751 0.9876013 -1.02820780
2 0.50    12571  1.625056 0.012361039 0.8740768  1.625751 0.9995725 -0.05622141
3 0.90    12730  1.621472 0.008907677 0.6246751  1.625751 0.9973680 -0.48037764
4 0.98    12836  1.628164 0.008420554 0.5810123  1.625751 1.0014844  0.28658810

The means agree to within 1.2 per cent across a walker that reverses at every step and one that travels almost in a straight line, and none sits more than 1.0 standard errors from the chord formula. The standard error comes from resampling animals, which is the safe default for a ratio of totals and which also allows for visits clustering within individuals. The coefficient of variation does not agree: it runs from 0.58 to 1.28. A tortuous animal delivers the same average visit through a wilder mixture of quick clips and long meanders, and that is variance you pay for in the interval, not in the estimate.

Where it earns its keep: resting animals

Now switch resting on. The animals are active 30 per cent of the time, in bouts of about a minute, with rests of about two minutes. Nothing else changes.

p_ra <- dt / 120                                  # mean rest bout 120 s
A_level <- 0.30
p_ar <- p_ra * (1 - A_level) / A_level            # gives 30 per cent active in the long run

set.seed(2703)
gr <- gas(nanim, nstep, dt, v, W, r, half, rho = 0.9, p_ar = p_ar, p_ra = p_ra)
lam_r <- length(gr$visits) / Tt
W_r <- mean(gr$visits)

rem_D <- function(y, t, vv, r, theta) y * pi / (t * vv * r * (2 + theta))
res <- data.frame(
  method = c("REST: (y/t) E[T] / a",
             "REM with the camera's active speed",
             "REM with the true time-averaged speed"),
  D_hat = c(lam_r * W_r / a_zone,
            rem_D(length(gr$visits), Tt, v, r, 2 * half),
            rem_D(length(gr$visits), Tt, v * mean(gr$active), r, 2 * half)))
res$ratio_to_truth <- res$D_hat / Dtrue
bs_r <- boot_stay(gr)
naive_se <- sd(gr$visits) / sqrt(length(gr$visits))    # pretends visits are independent
c(n_visits = length(gr$visits), realised_activity = mean(gr$active), mean_stay = W_r,
  boot_se = sd(bs_r), naive_se = naive_se, se_ratio = sd(bs_r) / naive_se,
  CV_stay = sd(gr$visits) / W_r, max_stay = max(gr$visits))

         n_visits realised_activity         mean_stay           boot_se 
     3114.0000000         0.3028683         5.2183915         0.5915992 
         naive_se          se_ratio           CV_stay          max_stay 
        0.5537578         1.0683357         5.9216465       695.7260233

res

                                 method      D_hat ratio_to_truth
1                  REST: (y/t) E[T] / a 0.29095640      0.9698547
2    REM with the camera's active speed 0.09064531      0.3021510
3 REM with the true time-averaged speed 0.29928950      0.9976317

This is the argument for the staying time model in one table. The animals now spend 70 per cent of their lives motionless, and the staying time model does not care: it returns 0.97 times the truth, with a bootstrap standard error of 0.11 on the ratio, without being told the activity level, because a resting animal in the zone is still an animal in the zone and Little’s law counts it. Feed the random encounter model the speed the camera actually measures, which is the speed of moving animals, and it returns 0.30 times the truth. It only recovers once you hand it the activity level from somewhere else.

The bill arrives in the shape of the distribution. The coefficient of variation of the staying time has gone from about 0.6 to 5.9, and the longest single visit was 696 seconds against a mean of 5.2. The mean of that distribution is a rare-event statistic: most visits are quick transits, and the average is set by the handful that turned into a nap. Even with 3114 visits recorded, the standard error on the mean is 11 per cent, and it carries straight through to the density.

One thing is not driving that: clustering. Resampling animals rather than visits changes the standard error by a factor of 1.07, so the visits are close to independent here and the imprecision is the tail itself, not the clustering. Between-camera heterogeneity is a different matter, and the checking post deals with it.

The camera has a stopwatch with a limit

Which brings the whole thing down. A camera records a clip of fixed length, or a burst with a maximum, and any stay that outlasts it is cut off. Cut a distribution whose mean lives in its tail and you cut the mean.

caps <- c(5, 10, 30, 60, 120, Inf)
cens <- data.frame(setting = ifelse(is.finite(caps), paste0(caps, " s"), "no limit"),
                   cap = caps)
cens$mean_recorded <- vapply(caps, function(cap) mean(pmin(gr$visits, cap)), numeric(1))
cens$D_hat <- lam_r * cens$mean_recorded / a_zone
cens$ratio_to_truth <- cens$D_hat / Dtrue
cens$pct_visits_cut <- vapply(caps, function(cap) 100 * mean(gr$visits > cap), numeric(1))
cens[, c("setting", "mean_recorded", "D_hat", "ratio_to_truth", "pct_visits_cut")]

   setting mean_recorded      D_hat ratio_to_truth pct_visits_cut
1      5 s      1.713229 0.09552274      0.3184091       3.050739
2     10 s      1.851698 0.10324320      0.3441440       2.697495
3     30 s      2.356079 0.13136543      0.4378848       2.312139
4     60 s      2.970055 0.16559824      0.5519941       1.830443
5    120 s      3.843034 0.21427203      0.7142401       1.027617
6 no limit      5.218391 0.29095640      0.9698547       0.000000

Line plot with maximum video length in seconds on the x axis, on a log scale from 5 to 120, and the estimated density on the y axis. The line rises steeply from a low value at 5 seconds and is still climbing towards the dashed horizontal truth line at 120 seconds. A second axis label notes the percentage of visits cut at each setting. — Figure 1: Density estimate from the staying time model against the maximum recording length, for a population that rests 70 per cent of the time. The dashed line is the truth. A camera set to a 10 second clip returns less than half the population, and the shortfall comes entirely from the small proportion of visits that get cut.

A 30 second clip recovers 44 per cent of the population while cutting only 2.3 per cent of the visits. That is the shape of the problem: the visits you lose are exactly the ones carrying the mean.

The instinct at this point is to reach for survival analysis, since these are textbook right-censored durations. It does not help, and the reason is worth understanding. When every observation is censored at the same known time $c$, the Kaplan-Meier curve is just the empirical survival function up to $c$ and stops. The restricted mean it gives you is $\mathbb{E}[\min(T, c)]$, which is exactly the truncated mean in the table above. The tail beyond $c$ is not poorly estimated, it is not estimated: nothing in the data speaks to it. You can bound $D$ from below, and that is all.

So the fix is not statistical, it is a camera setting, and the table tells you what to set it to. Extrapolating the tail with a parametric survival model is possible and is what a serious analysis does, but you are then reporting a density that depends on a distributional assumption about behaviour you did not observe, which should be said out loud rather than buried in a fitted object.

The honest ledger

The staying time model buys freedom from the speed measurement and pays for it in three places, all of which are real.

The zone area $a$ has to be measured, and it now sits in the denominator of the estimate alone rather than mixed into a perimeter. Rowcliffe et al. (2011) is the standard route to it, and it is not a tape measure job: the effective zone depends on the animal’s size and on the sensor, so it is a small distance sampling problem hiding inside the density estimate.

The staying time has to be captured whole, which the section above quantifies.

And the arrival has to be crisply defined, because $y$ counts entries into the zone and a camera counts triggers. Palencia et al. (2021) compare the staying time model against the alternatives on real data and find each has its own failure mode; the point of a simulation like this is not to crown a winner but to know which knob each one is sensitive to before you pick.

Where to go next

Both models now rest on quantities measured from the same cameras, in the same places, by the same field team. If those places are not a random sample of the landscape, everything above is precisely and confidently wrong. The checking post takes that apart, along with the interval, which is not what a Poisson count would suggest.

References

Nakashima, Fukasawa, Samejima 2018 Journal of Applied Ecology 55:735-744 (10.1111/1365-2664.13059)
Rowcliffe, Field, Turvey, Carbone 2008 Journal of Applied Ecology 45(4):1228-1236 (10.1111/j.1365-2664.2008.01473.x)
Hutchinson, Waser 2007 Biological Reviews 82(3):335-359 (10.1111/j.1469-185X.2007.00014.x)
Rowcliffe, Carbone, Jansen, Kays, Kranstauber 2011 Methods in Ecology and Evolution 2(5):464-476 (10.1111/j.2041-210X.2011.00094.x)
Palencia, Rowcliffe, Vicente, Acevedo 2021 Journal of Applied Ecology 58:1583-1592 (10.1111/1365-2664.13913)

--- title: "Staying time and the REST model in R" description: "The REST model estimates animal density from how long animals stay in view, with no external speed. Little's law and Cauchy's chord identity, tested in R." date: "2026-07-16 10:00" categories: [R, camera traps, abundance, movement ecology, ecology tutorial] image: thumbnail.png image-alt: "Plot of the density estimate against video length, showing the estimate falling steadily as the recording is cut shorter." --- The random encounter model has one weak joint: it needs the average speed, and the average has to include every hour the animal spends lying still. You cannot get that from the camera without an assumption about activity, and you cannot get it from a GPS collar without collaring something. The random encounter and staying time model (Nakashima et al. 2018) removes that joint. It uses the same cameras and the same triggers, but instead of asking how fast animals cross the detection zone it asks how long they stay in it. The estimator is: $$D = \frac{y}{t} \cdot \frac{\mathbb{E}[T_{\text{stay}}]}{a}$$ where $a$ is the area of the detection zone and $T_{\text{stay}}$ is the time an animal spends in view per visit. No speed. No activity level. Both are measured for you, implicitly, by the stopwatch. ## Why that works: Little's law The formula is not a camera trap result at all. It is Little's law, which says that for any stable system, the average number of customers inside equals the arrival rate times the average time each one spends there. Queueing theorists have used it since 1961 to reason about supermarket tills. Point it at a detection zone. The average number of animals inside is $D a$, because the gas is uniform. The arrival rate is $y/t$. The average time each one spends there is $\mathbb{E}[T_{\text{stay}}]$. So $$D a = \frac{y}{t}\,\mathbb{E}[T_{\text{stay}}]$$ and that is the whole derivation. Notice what is *not* in it: no assumption about how the animals move, no speed, no direction distribution, not even the shape of the zone. Little's law needs stationarity and nothing else. Set it beside the random encounter model, which says $y/t = D v r(2+\theta)/\pi$, and one substitution gives a prediction about staying time: $$\mathbb{E}[T_{\text{stay}}] = \frac{\pi a}{v L}, \qquad L = r(2+\theta)$$ That is Cauchy's mean chord formula wearing a wildlife hat: the average chord an isotropic straight line cuts across a convex region has length $\pi a / L$, and dividing by speed turns a length into a time. So the two models are the same equation read in two directions. The random encounter model measures $v$ elsewhere and predicts the count; the staying time model measures the count and the clock, and lets them tell it about $v$. ## The rig Same gas as the encounter rate post: independent walkers on a torus, exact boundary crossings, no time-step discretisation to apologise for. The extra piece here is the stopwatch. ```{r geometry} inside <- function(P, r, half) { (P[, 1]^2 + P[, 2]^2 <= r^2) & (abs(atan2(P[, 2], P[, 1])) <= half) } crossings <- function(P1, D, r, half) { m <- nrow(P1) U <- matrix(Inf, m, 4); IN <- matrix(FALSE, m, 4) a <- D[, 1]^2 + D[, 2]^2 b <- 2 * (P1[, 1] * D[, 1] + P1[, 2] * D[, 2]) cc <- P1[, 1]^2 + P1[, 2]^2 - r^2 disc <- b^2 - 4 * a * cc ok <- disc > 0 & a > 0 sq <- ifelse(ok, sqrt(pmax(disc, 0)), NA_real_) for (k in 1:2) { # the arc u <- (-b + c(-1, 1)[k] * sq) / (2 * a) Q <- cbind(P1[, 1] + u * D[, 1], P1[, 2] + u * D[, 2]) good <- ok & !is.na(u) & u > 0 & u <= 1 & abs(atan2(Q[, 2], Q[, 1])) <= half good[is.na(good)] <- FALSE U[good, k] <- u[good] IN[good, k] <- ((D[, 1] * Q[, 1] + D[, 2] * Q[, 2]) < 0)[good] } for (k in 1:2) { # the two radii s <- c(1, -1)[k] * half e <- c(cos(s), sin(s)) nrm <- c(cos(s + c(1, -1)[k] * pi / 2), sin(s + c(1, -1)[k] * pi / 2)) cD <- D[, 1] * e[2] - D[, 2] * e[1] u <- -(P1[, 1] * e[2] - P1[, 2] * e[1]) / cD tt <- (P1[, 1] + u * D[, 1]) * e[1] + (P1[, 2] + u * D[, 2]) * e[2] good <- abs(cD) > 1e-12 & !is.na(u) & u > 0 & u <= 1 & tt >= 0 & tt <= r good[is.na(good)] <- FALSE U[good, 2 + k] <- u[good] IN[good, 2 + k] <- ((D[, 1] * nrm[1] + D[, 2] * nrm[2]) < 0)[good] } list(U = U, IN = IN) } ``` An animal is active or resting, switching between the two as a Markov chain on the step grid. Resting means the step length is zero, which the crossing routine handles for free: a zero-length segment cannot cross anything. The proportion of time active follows from the two switching probabilities. ```{r walk} walk_one <- function(nstep, dt, v, W, r, half, rho, p_ar = 0, p_ra = 1) { ## active / resting state sequence, built bout by bout (p_ar = 0 means never rests) st <- if (p_ar <= 0) rep(1L, nstep) else { z <- integer(0); s <- 1L while (length(z) < nstep) { z <- c(z, rep(s, rgeom(1, if (s == 1L) p_ar else p_ra) + 1L)); s <- 1L - s } z[seq_len(nstep)] } u <- runif(nstep) turn <- 2 * atan(((1 - rho) / (1 + rho)) * tan(pi * (u - 0.5))) phi <- runif(1, 0, 2 * pi) + cumsum(turn) step <- v * dt * st # zero while resting D <- cbind(step * cos(phi), step * sin(phi)) p0 <- c(runif(1, -W / 2, W / 2), runif(1, -W / 2, W / 2)) X <- p0[1] + cumsum(D[, 1]); Y <- p0[2] + cumsum(D[, 2]) P1 <- cbind(c(p0[1], X[-nstep]), c(p0[2], Y[-nstep])) R1 <- cbind(((P1[, 1] + W / 2) %% W) - W / 2, ((P1[, 2] + W / 2) %% W) - W / 2) frac_in <- mean(inside(R1, r, half)) # time-average occupancy idx <- which((R1[, 1]^2 + R1[, 2]^2) < (r + v * dt)^2) if (!length(idx)) return(list(visits = numeric(0), frac_in = frac_in, active = mean(st))) cr <- crossings(R1[idx, , drop = FALSE], D[idx, , drop = FALSE], r, half) hit <- which(is.finite(cr$U), arr.ind = TRUE) if (!nrow(hit)) return(list(visits = numeric(0), frac_in = frac_in, active = mean(st))) tt <- (idx[hit[, 1]] - 1 + cr$U[hit]) * dt o <- order(tt); tt <- tt[o]; inw <- cr$IN[hit][o] t_in <- tt[inw]; t_out <- tt[!inw] if (!length(t_in)) return(list(visits = numeric(0), frac_in = frac_in, active = mean(st))) j <- findInterval(t_in, t_out) + 1L # the exit that closes each entry ok <- j <= length(t_out) # drop a visit still open at the end list(visits = t_out[j[ok]] - t_in[ok], frac_in = frac_in, active = mean(st)) } gas <- function(n, nstep, dt, v, W, r, half, rho, p_ar = 0, p_ra = 1) { out <- lapply(seq_len(n), function(i) walk_one(nstep, dt, v, W, r, half, rho, p_ar, p_ra)) list(visits = unlist(lapply(out, `[[`, "visits")), dur = vapply(out, function(z) sum(z$visits), numeric(1)), # per animal nvis = vapply(out, function(z) length(z$visits), numeric(1)), frac_in = vapply(out, `[[`, numeric(1), "frac_in"), active = vapply(out, `[[`, numeric(1), "active")) } ## the mean staying time is a ratio of two per-animal totals, so its standard error ## comes from resampling animals, which also allows for any clustering of visits boot_stay <- function(g, B = 600) { n <- length(g$dur) vapply(seq_len(B), function(b) { i <- sample.int(n, n, replace = TRUE) sum(g$dur[i]) / sum(g$nvis[i]) }, numeric(1)) } ``` ## Four identities, one simulation Start with a gas that never rests. There are four claims stacked on top of each other and the simulation checks each one separately: that the gas really is uniform, that Little's law holds, that the mean staying time matches Cauchy's chord formula, and that the resulting density estimate is right. ```{r pure} W <- 20; r <- 4; half <- 20 * pi / 180; v <- 1; dt <- 0.5 A <- W^2; a_zone <- half * r^2; L <- 2 * r + 2 * r * half nanim <- 120; nstep <- 20000; Tt <- nstep * dt; Dtrue <- nanim / A set.seed(2701) g <- gas(nanim, nstep, dt, v, W, r, half, rho = 0.9) lambda <- length(g$visits) / Tt # arrivals per second, whole population Wbar <- mean(g$visits) # mean staying time Lnum <- sum(g$frac_in) # mean number inside, from positions data.frame( claim = c("gas is uniform: mean number inside = D a", "Little's law: arrival rate x staying time = mean number inside", "Cauchy chord: mean staying time = pi a / (v L)", "REST estimate of density"), measured = c(Lnum, lambda * Wbar, Wbar, lambda * Wbar / a_zone), predicted = c(Dtrue * a_zone, Lnum, pi * a_zone / (v * L), Dtrue)) ``` Every row lands. The third is the one worth pausing on: the mean staying time is `r sprintf("%.3f", Wbar)` seconds against a predicted `r sprintf("%.3f", pi * a_zone / (v * L))`, and that prediction came from a nineteenth-century integral geometry result about random lines crossing a convex set. Now vary the tortuosity. The encounter rate post showed that path persistence does not change how often animals arrive. Little's law then forces something less obvious: it cannot change how long they stay either, because the occupancy and the arrival rate are both fixed, and the staying time is their ratio. ```{r tort} set.seed(2702) tor <- do.call(rbind, lapply(c(0, 0.5, 0.9, 0.98), function(rh) { gg <- gas(150, nstep, dt, v, W, r, half, rho = rh) bs <- boot_stay(gg) data.frame(rho = rh, n_visits = length(gg$visits), mean_stay = mean(gg$visits), boot_se = sd(bs), CV_stay = sd(gg$visits) / mean(gg$visits)) })) tor$predicted <- pi * a_zone / (v * L) tor$ratio <- tor$mean_stay / tor$predicted tor$z <- (tor$mean_stay - tor$predicted) / tor$boot_se tor ``` The means agree to within `r sprintf("%.1f", 100 * max(abs(tor$ratio - 1)))` per cent across a walker that reverses at every step and one that travels almost in a straight line, and none sits more than `r sprintf("%.1f", max(abs(tor$z)))` standard errors from the chord formula. The standard error comes from resampling animals, which is the safe default for a ratio of totals and which also allows for visits clustering within individuals. The coefficient of variation does not agree: it runs from `r sprintf("%.2f", min(tor$CV_stay))` to `r sprintf("%.2f", max(tor$CV_stay))`. A tortuous animal delivers the same average visit through a wilder mixture of quick clips and long meanders, and that is variance you pay for in the interval, not in the estimate. ## Where it earns its keep: resting animals Now switch resting on. The animals are active `r sprintf("%.0f", 30)` per cent of the time, in bouts of about a minute, with rests of about two minutes. Nothing else changes. ```{r resting} p_ra <- dt / 120 # mean rest bout 120 s A_level <- 0.30 p_ar <- p_ra * (1 - A_level) / A_level # gives 30 per cent active in the long run set.seed(2703) gr <- gas(nanim, nstep, dt, v, W, r, half, rho = 0.9, p_ar = p_ar, p_ra = p_ra) lam_r <- length(gr$visits) / Tt W_r <- mean(gr$visits) rem_D <- function(y, t, vv, r, theta) y * pi / (t * vv * r * (2 + theta)) res <- data.frame( method = c("REST: (y/t) E[T] / a", "REM with the camera's active speed", "REM with the true time-averaged speed"), D_hat = c(lam_r * W_r / a_zone, rem_D(length(gr$visits), Tt, v, r, 2 * half), rem_D(length(gr$visits), Tt, v * mean(gr$active), r, 2 * half))) res$ratio_to_truth <- res$D_hat / Dtrue bs_r <- boot_stay(gr) naive_se <- sd(gr$visits) / sqrt(length(gr$visits)) # pretends visits are independent c(n_visits = length(gr$visits), realised_activity = mean(gr$active), mean_stay = W_r, boot_se = sd(bs_r), naive_se = naive_se, se_ratio = sd(bs_r) / naive_se, CV_stay = sd(gr$visits) / W_r, max_stay = max(gr$visits)) res ``` This is the argument for the staying time model in one table. The animals now spend `r sprintf("%.0f", 100 * (1 - mean(gr$active)))` per cent of their lives motionless, and the staying time model does not care: it returns `r sprintf("%.2f", res$ratio_to_truth[1])` times the truth, with a bootstrap standard error of `r sprintf("%.2f", sd(bs_r) / W_r)` on the ratio, without being told the activity level, because a resting animal in the zone is still an animal in the zone and Little's law counts it. Feed the random encounter model the speed the camera actually measures, which is the speed of moving animals, and it returns `r sprintf("%.2f", res$ratio_to_truth[2])` times the truth. It only recovers once you hand it the activity level from somewhere else. The bill arrives in the shape of the distribution. The coefficient of variation of the staying time has gone from about `r sprintf("%.1f", tor$CV_stay[tor$rho == 0.9])` to `r sprintf("%.1f", sd(gr$visits) / W_r)`, and the longest single visit was `r sprintf("%.0f", max(gr$visits))` seconds against a mean of `r sprintf("%.1f", W_r)`. The mean of that distribution is a rare-event statistic: most visits are quick transits, and the average is set by the handful that turned into a nap. Even with `r sprintf("%d", length(gr$visits))` visits recorded, the standard error on the mean is `r sprintf("%.0f", 100 * sd(bs_r) / W_r)` per cent, and it carries straight through to the density. One thing is not driving that: clustering. Resampling animals rather than visits changes the standard error by a factor of `r sprintf("%.2f", sd(bs_r) / naive_se)`, so the visits are close to independent here and the imprecision is the tail itself, not the clustering. Between-camera heterogeneity is a different matter, and the checking post deals with it. ## The camera has a stopwatch with a limit Which brings the whole thing down. A camera records a clip of fixed length, or a burst with a maximum, and any stay that outlasts it is cut off. Cut a distribution whose mean lives in its tail and you cut the mean. ```{r censor} caps <- c(5, 10, 30, 60, 120, Inf) cens <- data.frame(setting = ifelse(is.finite(caps), paste0(caps, " s"), "no limit"), cap = caps) cens$mean_recorded <- vapply(caps, function(cap) mean(pmin(gr$visits, cap)), numeric(1)) cens$D_hat <- lam_r * cens$mean_recorded / a_zone cens$ratio_to_truth <- cens$D_hat / Dtrue cens$pct_visits_cut <- vapply(caps, function(cap) 100 * mean(gr$visits > cap), numeric(1)) cens[, c("setting", "mean_recorded", "D_hat", "ratio_to_truth", "pct_visits_cut")] ``` ```{r fig-censoring} #| echo: false #| fig-cap: "Density estimate from the staying time model against the maximum recording length, for a population that rests 70 per cent of the time. The dashed line is the truth. A camera set to a 10 second clip returns less than half the population, and the shortfall comes entirely from the small proportion of visits that get cut." #| fig-alt: "Line plot with maximum video length in seconds on the x axis, on a log scale from 5 to 120, and the estimated density on the y axis. The line rises steeply from a low value at 5 seconds and is still climbing towards the dashed horizontal truth line at 120 seconds. A second axis label notes the percentage of visits cut at each setting." #| warning: false library(ggplot2) fc <- cens[is.finite(cens$cap), ] ggplot(fc, aes(cap, D_hat)) + geom_hline(yintercept = Dtrue, linetype = 2, colour = "#93a87f") + geom_line(colour = "#275139", linewidth = 0.8) + geom_point(size = 2.6, colour = "#275139") + geom_text(aes(label = sprintf("%.1f%% cut", pct_visits_cut)), vjust = 1.9, size = 3.1, colour = "#5d6b61") + annotate("text", x = 6, y = Dtrue * 1.03, hjust = 0, size = 3.3, colour = "#46604a", label = "true density") + scale_x_log10(breaks = fc$cap) + coord_cartesian(ylim = c(0, Dtrue * 1.12)) + labs(x = "maximum recording length (s)", y = "estimated density (per m2)") + theme_minimal(base_size = 12) + theme(panel.grid.minor = element_blank()) ``` A 30 second clip recovers `r sprintf("%.0f", 100 * cens$ratio_to_truth[cens$cap == 30])` per cent of the population while cutting only `r sprintf("%.1f", cens$pct_visits_cut[cens$cap == 30])` per cent of the visits. That is the shape of the problem: the visits you lose are exactly the ones carrying the mean. The instinct at this point is to reach for survival analysis, since these are textbook right-censored durations. It does not help, and the reason is worth understanding. When every observation is censored at the same known time $c$, the Kaplan-Meier curve is just the empirical survival function up to $c$ and stops. The restricted mean it gives you is $\mathbb{E}[\min(T, c)]$, which is exactly the truncated mean in the table above. The tail beyond $c$ is not poorly estimated, it is not estimated: nothing in the data speaks to it. You can bound $D$ from below, and that is all. So the fix is not statistical, it is a camera setting, and the table tells you what to set it to. Extrapolating the tail with a parametric survival model is possible and is what a serious analysis does, but you are then reporting a density that depends on a distributional assumption about behaviour you did not observe, which should be said out loud rather than buried in a fitted object. ## The honest ledger The staying time model buys freedom from the speed measurement and pays for it in three places, all of which are real. The zone area $a$ has to be measured, and it now sits in the denominator of the estimate alone rather than mixed into a perimeter. Rowcliffe et al. (2011) is the standard route to it, and it is not a tape measure job: the effective zone depends on the animal's size and on the sensor, so it is a small distance sampling problem hiding inside the density estimate. The staying time has to be captured whole, which the section above quantifies. And the arrival has to be crisply defined, because $y$ counts entries into the zone and a camera counts triggers. Palencia et al. (2021) compare the staying time model against the alternatives on real data and find each has its own failure mode; the point of a simulation like this is not to crown a winner but to know which knob each one is sensitive to before you pick. ## Where to go next Both models now rest on quantities measured from the same cameras, in the same places, by the same field team. If those places are not a random sample of the landscape, everything above is precisely and confidently wrong. The checking post takes that apart, along with the interval, which is not what a Poisson count would suggest. ## Related tutorials - [Camera trap density: the random encounter model](../random-encounter-model-camera-traps/) - [Measuring animal speed from camera traps](../animal-speed-from-camera-traps/) - [Kaplan-Meier survival curves in R](../kaplan-meier-survival-curves/) - [Checking a camera trap density estimate](../checking-a-camera-trap-density-estimate/) ## References - Nakashima, Fukasawa, Samejima 2018 Journal of Applied Ecology 55:735-744 (10.1111/1365-2664.13059) - Rowcliffe, Field, Turvey, Carbone 2008 Journal of Applied Ecology 45(4):1228-1236 (10.1111/j.1365-2664.2008.01473.x) - Hutchinson, Waser 2007 Biological Reviews 82(3):335-359 (10.1111/j.1469-185X.2007.00014.x) - Rowcliffe, Carbone, Jansen, Kays, Kranstauber 2011 Methods in Ecology and Evolution 2(5):464-476 (10.1111/j.2041-210X.2011.00094.x) - Palencia, Rowcliffe, Vicente, Acevedo 2021 Journal of Applied Ecology 58:1583-1592 (10.1111/1365-2664.13913)