Step selection functions for animal movement

movement ecology

habitat selection

ecology tutorial

Fit step selection functions in R from scratch: matched steps, a conditional logistic likelihood by hand, and why pooling strata attenuates selection.

Author

Tidy Ecology

Published

2026-07-13

A resource selection function asks whether the places an animal used hold more of a resource than the places available to it, with availability sampled across a home range or study area. That framing hides a problem: an animal cannot reach every cell in one step. From where it stands now, only nearby cells are within a single move, so “available” should really mean “reachable in one step”. Step selection functions make availability local. For each observed step, they compare the resource at the step the animal took against the resource at a set of steps it could have taken from the same starting point. The machinery that makes this comparison honest is a conditional (stratified) logistic likelihood, and we build it here from nothing but optim.

A track that selects for a resource

We work on a smooth resource surface and simulate a track whose owner genuinely prefers rich cells. At each step the animal draws a pool of candidate moves from a movement kernel (a gamma step length and a von Mises turning angle relative to its current heading), then picks one candidate with probability proportional to exp(beta * resource) at its endpoint. This is the step selection process written as a generative model, so the true selection strength is a known number.

library(ggplot2)
library(grid)

res_raw <- function(x, y) {
  1.5 * exp(-((x - 38)^2 + (y - 62)^2) / (2 * 28^2)) +
  1.2 * exp(-((x - 70)^2 + (y - 40)^2) / (2 * 32^2)) +
  (-0.00030) * ((x - 50)^2 + (y - 50)^2)
}
gx <- seq(0, 100, length.out = 120)
grid_df <- expand.grid(x = gx, y = gx)
grid_df$raw <- res_raw(grid_df$x, grid_df$y)
mu_r <- mean(grid_df$raw); sd_r <- sd(grid_df$raw)
resource <- function(x, y) (res_raw(x, y) - mu_r) / sd_r
grid_df$z <- resource(grid_df$x, grid_df$y)

We need to draw turning angles, so we write a von Mises sampler by hand (the Best and Fisher 1979 rejection method). The same distribution supplied the emissions in the movement HMM; here it shapes where the available steps point.

rvm <- function(n, mu, kappa) {
  out <- numeric(n)
  a <- 1 + sqrt(1 + 4 * kappa^2)
  b <- (a - sqrt(2 * a)) / (2 * kappa)
  r <- (1 + b^2) / (2 * b)
  for (i in seq_len(n)) {
    repeat {
      u1 <- runif(1); u2 <- runif(1); u3 <- runif(1)
      z <- cos(pi * u1)
      f <- (1 + r * z) / (r + z)
      cc <- kappa * (r - f)
      if (cc * (2 - cc) - u2 > 0 || log(cc / u2) + 1 - cc >= 0) break
    }
    theta <- sign(u3 - 0.5) * acos(f)
    out[i] <- theta + mu
  }
  ((out + pi) %% (2 * pi)) - pi
}

sh_move <- 2.2; sc_move <- 5.5; kap_move <- 0.9   # movement kernel (truth)
beta_true <- 0.85                                  # selection strength (truth)
n <- 1200; M <- 300                                # steps; candidate pool per step

set.seed(4711)
pos <- matrix(NA_real_, n + 1, 2); pos[1, ] <- c(45, 50)
heading <- numeric(n + 1); heading[1] <- runif(1, -pi, pi)
for (i in 1:n) {
  L <- rgamma(M, shape = sh_move, scale = sc_move)
  a <- rvm(M, 0, kap_move)
  ang <- heading[i] + a
  ex <- pos[i, 1] + L * cos(ang)
  ey <- pos[i, 2] + L * sin(ang)
  w <- exp(beta_true * resource(ex, ey))
  j <- sample.int(M, 1, prob = w)
  pos[i + 1, ] <- c(ex[j], ey[j]); heading[i + 1] <- ang[j]
}
trk <- data.frame(x = pos[, 1], y = pos[, 2])

dx <- diff(pos[, 1]); dy <- diff(pos[, 2])
step_len <- sqrt(dx^2 + dy^2)
bearing <- atan2(dy, dx)
turn <- ((diff(bearing) + pi) %% (2 * pi)) - pi
mean_step <- mean(step_len); med_step <- median(step_len)

The track has 1200 steps with a mean length of 11.5 units (median 9.9). To sample available steps we need a tentative movement kernel, so we fit a gamma to the observed step lengths by matching moments and a von Mises to the observed turns through the mean resultant length.

m1 <- mean(step_len); v1 <- var(step_len)
sh_hat <- m1^2 / v1; sc_hat <- v1 / m1
C <- mean(cos(turn)); S <- mean(sin(turn)); Rbar <- sqrt(C^2 + S^2)
mu_turn <- atan2(S, C)
kap_hat <- if (Rbar < 0.53) 2 * Rbar + Rbar^3 + 5 * Rbar^5 / 6 else
           if (Rbar < 0.85) -0.4 + 1.39 * Rbar + 0.43 / (1 - Rbar) else
           1 / (Rbar^3 - 4 * Rbar^2 + 3 * Rbar)

Matched steps: one used, ten available

Now the core construction. Each observed step becomes a stratum. The used step (the move the animal made) is the case; alongside it we place K available steps that start from the same point, with lengths from the tentative gamma and turning angles from the tentative von Mises. Every stratum shares its starting location, so the resource that happens to sit near that spot is common to the used and the available steps in that stratum.

K <- 10
set.seed(20260713)
use_idx <- 2:n
n_str <- length(use_idx)
Zcase <- resource(pos[use_idx + 1, 1], pos[use_idx + 1, 2])
Zctrl <- matrix(NA_real_, n_str, K)
for (s in seq_along(use_idx)) {
  i <- use_idx[s]
  Ls <- rgamma(K, shape = sh_hat, scale = sc_hat)
  as <- rvm(K, mu_turn, kap_hat)
  ang <- heading[i] + as
  Zctrl[s, ] <- resource(pos[i, 1] + Ls * cos(ang), pos[i, 2] + Ls * sin(ang))
}
Z <- cbind(Zcase, Zctrl)          # n_str x (K+1); column 1 is the used step

d2 <- (pos[use_idx, 1] - 40)^2 + (pos[use_idx, 2] - 60)^2
ex_s <- which.min(d2); i_ex <- use_idx[ex_s]
set.seed(77)
Le <- rgamma(K, shape = sh_hat, scale = sc_hat)
ae <- rvm(K, mu_turn, kap_hat)
ange <- heading[i_ex] + ae
fan <- data.frame(x0 = pos[i_ex, 1], y0 = pos[i_ex, 2],
                  x1 = pos[i_ex, 1] + Le * cos(ange), y1 = pos[i_ex, 2] + Le * sin(ange))
used_seg <- data.frame(x0 = pos[i_ex, 1], y0 = pos[i_ex, 2],
                       x1 = pos[i_ex + 1, 1], y1 = pos[i_ex + 1, 2])
ggplot() +
  geom_raster(data = grid_df, aes(x, y, fill = z)) +
  scale_fill_gradient(low = "#f5f4ee", high = "#275139", name = "resource\n(SD)") +
  geom_path(data = trk, aes(x, y), colour = "#16241d", linewidth = 0.3, alpha = 0.7) +
  geom_segment(data = fan, aes(x0, y0, xend = x1, yend = y1),
               colour = "#5d6b61", linewidth = 0.4, alpha = 0.9) +
  geom_point(data = fan, aes(x1, y1), colour = "#5d6b61", size = 1.3) +
  geom_segment(data = used_seg, aes(x0, y0, xend = x1, yend = y1), colour = "#b5534e",
               linewidth = 0.9, arrow = arrow(length = unit(0.16, "cm"), type = "closed")) +
  geom_point(data = fan, aes(x0, y0), colour = "#16241d", size = 1.8) +
  coord_equal(xlim = c(0, 100), ylim = c(0, 100), expand = FALSE) +
  labs(x = "easting", y = "northing") +
  theme_minimal(base_size = 11) +
  theme(panel.grid = element_blank(), axis.text = element_text(colour = "#46604a"),
        legend.title = element_text(size = 8, colour = "#46604a"))

A green-shaded resource surface with a thin dark track winding across it. Near a bright patch, a red arrow marks the step the animal took and ten grey line segments show alternative steps radiating from the same point. — Figure 1: A simulated track on a resource surface (paper to forest is low to high). One step is shown as the used move (red arrow) against ten available moves (grey) that fan out from the same starting point. Each such set is one stratum in the conditional model.

The conditional likelihood, by hand

Within a stratum the animal chose the used step out of the K + 1 on offer. If selection is proportional to exp(z * beta), the probability that this particular step was the one taken is a softmax over the stratum:

\[ \Pr(\text{used}_s) = \frac{\exp(z_{s,\text{used}}\,\beta)}{\sum_{j} \exp(z_{s,j}\,\beta)} . \]

The starting location enters the numerator and every term of the denominator identically, so it cancels. That is the whole point: the local abundance of the resource, which differs from stratum to stratum, drops out, and what remains is the selection for the resource over and above what was locally available. Summing the log of this over strata gives the conditional log-likelihood, which we hand to optim.

clogit_nll <- function(b) {
  eta <- Z * b
  -sum(eta[, 1] - log(rowSums(exp(eta))))
}
opt <- optim(0, clogit_nll, method = "BFGS", hessian = TRUE,
             control = list(reltol = 1e-10))
b_hat <- opt$par
se_hat <- sqrt(1 / opt$hessian[1, 1])
ci <- b_hat + c(-1, 1) * 1.96 * se_hat
rss <- exp(b_hat)

The estimate is 0.741 with a standard error of 0.059 and a 95% interval of [0.63, 0.86]. Exponentiating gives the relative selection strength: exp(beta) is 2.1, so a step ending in a cell one standard deviation richer in the resource is about 2.1 times as likely as an otherwise identical step ending in an average cell. This estimator is a conditional logistic regression; equivalently it is a Cox proportional hazards model stratified by step, which is how most software fits it.

Why the strata are not optional

It is tempting to throw the used and available steps into one ordinary logistic regression with a single intercept and forget the pairing. That intercept has to stand for the average availability across the whole track, but availability is not constant: as the animal moves through richer and poorer parts of the surface, the resource near its current position rises and falls. A single intercept cannot follow it, and the selection coefficient absorbs the damage.

long_z <- as.vector(t(Z))
long_y <- rep(c(1, rep(0, K)), n_str)
g_pool <- glm(long_y ~ long_z, family = binomial)
b_pool <- unname(coef(g_pool)[2]); se_pool <- summary(g_pool)$coefficients[2, 2]

Pooling collapses the coefficient to 0.142, roughly a fifth of the conditional estimate and close to no selection at all. The signal has not weakened; it has been averaged away by an intercept that could not track the moving baseline. Conditioning on the stratum is what protects it.

Recovering the truth, and where the bias comes from

Because this is a simulation we know the answer, so we can check the method directly. If we draw the available steps from the movement kernel that actually generated the track, rather than from a kernel fitted to the observed steps, the estimate should land on the true value.

set.seed(424242)
Zctrl_true <- matrix(NA_real_, n_str, K)
for (s in seq_along(use_idx)) {
  i <- use_idx[s]
  Ls <- rgamma(K, shape = sh_move, scale = sc_move)
  as <- rvm(K, 0, kap_move)
  ang <- heading[i] + as
  Zctrl_true[s, ] <- resource(pos[i, 1] + Ls * cos(ang), pos[i, 2] + Ls * sin(ang))
}
Zt <- cbind(Zcase, Zctrl_true)
clogit_nll_t <- function(b) { eta <- Zt * b; -sum(eta[, 1] - log(rowSums(exp(eta)))) }
opt_t <- optim(0, clogit_nll_t, method = "BFGS", hessian = TRUE, control = list(reltol = 1e-10))
b_hat_true <- opt_t$par; se_hat_true <- sqrt(1 / opt_t$hessian[1, 1])

With the true movement kernel supplying availability the estimate is 0.854, sitting on the generative value of 0.85. The fitted-kernel estimate above (0.74) is a shade lower, and the reason is instructive: the tentative gamma and von Mises were fitted to the observed steps, but the observed steps are already the product of selection, so the availability distribution they imply is slightly wrong. Estimating the movement kernel and the habitat selection together, rather than fixing one and inferring the other, removes this bias. That is integrated step selection analysis, and it is the next post.

dens_df <- rbind(data.frame(kind = "used steps", z = Zcase),
                 data.frame(kind = "available steps", z = as.vector(Zctrl)))
pa <- ggplot(dens_df, aes(z, fill = kind, colour = kind)) +
  geom_density(alpha = 0.45, linewidth = 0.5) +
  scale_fill_manual(values = c("used steps" = "#275139", "available steps" = "#93a87f")) +
  scale_colour_manual(values = c("used steps" = "#275139", "available steps" = "#93a87f")) +
  labs(x = "resource at step endpoint (SD)", y = "density", fill = NULL, colour = NULL) +
  theme_minimal(base_size = 11) +
  theme(panel.grid.minor = element_blank(), legend.position = c(0.72, 0.85),
        legend.background = element_rect(fill = "#ffffff", colour = NA))
coef_df <- data.frame(
  model = factor(c("conditional: movement kernel", "conditional: fitted kernel",
                   "pooled: strata ignored"),
                 levels = c("pooled: strata ignored", "conditional: fitted kernel",
                            "conditional: movement kernel")),
  est = c(b_hat_true, b_hat, b_pool), se = c(se_hat_true, se_hat, se_pool))
pb <- ggplot(coef_df, aes(est, model)) +
  geom_vline(xintercept = beta_true, linetype = "dashed", colour = "#b5534e") +
  geom_errorbarh(aes(xmin = est - 1.96 * se, xmax = est + 1.96 * se),
                 height = 0.16, colour = "#275139", linewidth = 0.6) +
  geom_point(colour = "#275139", size = 2.6) +
  annotate("text", x = beta_true, y = 3.44, label = "true beta",
           colour = "#b5534e", size = 3, hjust = -0.1) +
  labs(x = "selection coefficient", y = NULL) +
  theme_minimal(base_size = 11) + theme(panel.grid.minor = element_blank())

Warning: `geom_errorbarh()` was deprecated in ggplot2 4.0.0.
ℹ Please use the `orientation` argument of `geom_errorbar()` instead.

grid.newpage()
pushViewport(viewport(layout = grid.layout(1, 2, widths = unit(c(1, 1), "null"))))
print(pa, vp = viewport(layout.pos.row = 1, layout.pos.col = 1))
print(pb, vp = viewport(layout.pos.row = 1, layout.pos.col = 2))

`height` was translated to `width`.

Two panels. On the left, two overlapping density curves: used steps shifted right of available steps on the resource axis. On the right, three coefficient estimates with error bars against a dashed vertical line at the true value; the two conditional estimates sit near the line while the pooled estimate sits far to the left near zero. — Figure 2: Left: the resource at used steps sits higher than at available steps, the selection signal. Right: the conditional estimate recovers the true coefficient (dashed) when availability comes from the movement kernel, is slightly low with a fitted kernel, and collapses toward zero when the strata are pooled into one logistic regression.

Selection at which scale

A step selection function and a resource selection function answer different questions, and the numbers show it. Fitting a plain resource selection function here, with available points drawn from across the whole domain rather than from reachable steps, gives a coefficient of a different size again.

set.seed(999)
used_pts <- trk[-1, ]
n_av <- nrow(used_pts) * K
av_pts <- data.frame(x = runif(n_av, 0, 100), y = runif(n_av, 0, 100))
rsf_y <- c(rep(1, nrow(used_pts)), rep(0, n_av))
rsf_z <- c(resource(used_pts$x, used_pts$y), resource(av_pts$x, av_pts$y))
g_rsf <- glm(rsf_y ~ rsf_z, family = binomial)
b_rsf <- unname(coef(g_rsf)[2])

The domain-wide resource selection function returns 0.65. It is not wrong, and it is not the same thing: it measures broad-scale use of the landscape, the tendency to end up in rich areas over the whole track, whereas the step selection coefficient measures the choice made at each move given where the animal already was. These are different orders of selection in the sense of Johnson (1980), and a study should be clear about which one it is asking for. The step scale is where movement and habitat interact, and it is where the next post puts them together.

References

Fortin D, Beyer HL, Boyce MS, Smith DW, Duchesne T, Mao JS 2005. Ecology 86(5):1320-1330 (10.1890/04-0953).

Thurfjell H, Ciuti S, Boyce MS 2014. Movement Ecology 2:4 (10.1186/s40462-014-0004-7).

Forester JD, Im HK, Rathouz PJ 2009. Ecology 90(12):3554-3565 (10.1890/08-0874.1).

Johnson DH 1980. Ecology 61(1):65-71 (10.2307/1937156).

Manly BFJ, McDonald LL, Thomas DL, McDonald TL, Erickson WP 2002. Resource Selection by Animals, 2nd edn. Kluwer (ISBN 978-1-4020-0677-7).

Best DJ, Fisher NI 1979. Journal of the Royal Statistical Society Series C 28(2):152-157 (10.2307/2346732).

--- title: "Step selection functions for animal movement" description: "Fit step selection functions in R from scratch: matched steps, a conditional logistic likelihood by hand, and why pooling strata attenuates selection." date: "2026-07-13 09:00" categories: [movement ecology, habitat selection, R, ecology tutorial] image: thumbnail.png image-alt: "A resource surface with a simulated animal track and one step shown as a used move against a fan of available moves from the same start." --- A [resource selection function](../resource-selection-functions/) asks whether the places an animal used hold more of a resource than the places available to it, with availability sampled across a home range or study area. That framing hides a problem: an animal cannot reach every cell in one step. From where it stands now, only nearby cells are within a single move, so "available" should really mean "reachable in one step". Step selection functions make availability local. For each observed step, they compare the resource at the step the animal took against the resource at a set of steps it could have taken from the same starting point. The machinery that makes this comparison honest is a conditional (stratified) logistic likelihood, and we build it here from nothing but `optim`. ## A track that selects for a resource We work on a smooth resource surface and simulate a track whose owner genuinely prefers rich cells. At each step the animal draws a pool of candidate moves from a movement kernel (a gamma step length and a von Mises turning angle relative to its current heading), then picks one candidate with probability proportional to `exp(beta * resource)` at its endpoint. This is the step selection process written as a generative model, so the true selection strength is a known number. ```{r} #| label: setup #| message: false #| warning: false library(ggplot2) library(grid) res_raw <- function(x, y) { 1.5 * exp(-((x - 38)^2 + (y - 62)^2) / (2 * 28^2)) + 1.2 * exp(-((x - 70)^2 + (y - 40)^2) / (2 * 32^2)) + (-0.00030) * ((x - 50)^2 + (y - 50)^2) } gx <- seq(0, 100, length.out = 120) grid_df <- expand.grid(x = gx, y = gx) grid_df$raw <- res_raw(grid_df$x, grid_df$y) mu_r <- mean(grid_df$raw); sd_r <- sd(grid_df$raw) resource <- function(x, y) (res_raw(x, y) - mu_r) / sd_r grid_df$z <- resource(grid_df$x, grid_df$y) ``` We need to draw turning angles, so we write a von Mises sampler by hand (the Best and Fisher 1979 rejection method). The same distribution supplied the emissions in the [movement HMM](../two-state-movement-hmm/); here it shapes where the available steps point. ```{r} #| label: vonmises rvm <- function(n, mu, kappa) { out <- numeric(n) a <- 1 + sqrt(1 + 4 * kappa^2) b <- (a - sqrt(2 * a)) / (2 * kappa) r <- (1 + b^2) / (2 * b) for (i in seq_len(n)) { repeat { u1 <- runif(1); u2 <- runif(1); u3 <- runif(1) z <- cos(pi * u1) f <- (1 + r * z) / (r + z) cc <- kappa * (r - f) if (cc * (2 - cc) - u2 > 0 || log(cc / u2) + 1 - cc >= 0) break } theta <- sign(u3 - 0.5) * acos(f) out[i] <- theta + mu } ((out + pi) %% (2 * pi)) - pi } ``` ```{r} #| label: simulate sh_move <- 2.2; sc_move <- 5.5; kap_move <- 0.9 # movement kernel (truth) beta_true <- 0.85 # selection strength (truth) n <- 1200; M <- 300 # steps; candidate pool per step set.seed(4711) pos <- matrix(NA_real_, n + 1, 2); pos[1, ] <- c(45, 50) heading <- numeric(n + 1); heading[1] <- runif(1, -pi, pi) for (i in 1:n) { L <- rgamma(M, shape = sh_move, scale = sc_move) a <- rvm(M, 0, kap_move) ang <- heading[i] + a ex <- pos[i, 1] + L * cos(ang) ey <- pos[i, 2] + L * sin(ang) w <- exp(beta_true * resource(ex, ey)) j <- sample.int(M, 1, prob = w) pos[i + 1, ] <- c(ex[j], ey[j]); heading[i + 1] <- ang[j] } trk <- data.frame(x = pos[, 1], y = pos[, 2]) dx <- diff(pos[, 1]); dy <- diff(pos[, 2]) step_len <- sqrt(dx^2 + dy^2) bearing <- atan2(dy, dx) turn <- ((diff(bearing) + pi) %% (2 * pi)) - pi mean_step <- mean(step_len); med_step <- median(step_len) ``` The track has `r n` steps with a mean length of `r round(mean_step, 1)` units (median `r round(med_step, 1)`). To sample available steps we need a tentative movement kernel, so we fit a gamma to the observed step lengths by matching moments and a von Mises to the observed turns through the mean resultant length. ```{r} #| label: tentative m1 <- mean(step_len); v1 <- var(step_len) sh_hat <- m1^2 / v1; sc_hat <- v1 / m1 C <- mean(cos(turn)); S <- mean(sin(turn)); Rbar <- sqrt(C^2 + S^2) mu_turn <- atan2(S, C) kap_hat <- if (Rbar < 0.53) 2 * Rbar + Rbar^3 + 5 * Rbar^5 / 6 else if (Rbar < 0.85) -0.4 + 1.39 * Rbar + 0.43 / (1 - Rbar) else 1 / (Rbar^3 - 4 * Rbar^2 + 3 * Rbar) ``` ## Matched steps: one used, ten available Now the core construction. Each observed step becomes a stratum. The used step (the move the animal made) is the case; alongside it we place `K` available steps that start from the same point, with lengths from the tentative gamma and turning angles from the tentative von Mises. Every stratum shares its starting location, so the resource that happens to sit near that spot is common to the used and the available steps in that stratum. ```{r} #| label: strata K <- 10 set.seed(20260713) use_idx <- 2:n n_str <- length(use_idx) Zcase <- resource(pos[use_idx + 1, 1], pos[use_idx + 1, 2]) Zctrl <- matrix(NA_real_, n_str, K) for (s in seq_along(use_idx)) { i <- use_idx[s] Ls <- rgamma(K, shape = sh_hat, scale = sc_hat) as <- rvm(K, mu_turn, kap_hat) ang <- heading[i] + as Zctrl[s, ] <- resource(pos[i, 1] + Ls * cos(ang), pos[i, 2] + Ls * sin(ang)) } Z <- cbind(Zcase, Zctrl) # n_str x (K+1); column 1 is the used step ``` ```{r} #| label: fig-strata #| fig-cap: "A simulated track on a resource surface (paper to forest is low to high). One step is shown as the used move (red arrow) against ten available moves (grey) that fan out from the same starting point. Each such set is one stratum in the conditional model." #| fig-alt: "A green-shaded resource surface with a thin dark track winding across it. Near a bright patch, a red arrow marks the step the animal took and ten grey line segments show alternative steps radiating from the same point." #| fig-width: 6.4 #| fig-height: 5.0 d2 <- (pos[use_idx, 1] - 40)^2 + (pos[use_idx, 2] - 60)^2 ex_s <- which.min(d2); i_ex <- use_idx[ex_s] set.seed(77) Le <- rgamma(K, shape = sh_hat, scale = sc_hat) ae <- rvm(K, mu_turn, kap_hat) ange <- heading[i_ex] + ae fan <- data.frame(x0 = pos[i_ex, 1], y0 = pos[i_ex, 2], x1 = pos[i_ex, 1] + Le * cos(ange), y1 = pos[i_ex, 2] + Le * sin(ange)) used_seg <- data.frame(x0 = pos[i_ex, 1], y0 = pos[i_ex, 2], x1 = pos[i_ex + 1, 1], y1 = pos[i_ex + 1, 2]) ggplot() + geom_raster(data = grid_df, aes(x, y, fill = z)) + scale_fill_gradient(low = "#f5f4ee", high = "#275139", name = "resource\n(SD)") + geom_path(data = trk, aes(x, y), colour = "#16241d", linewidth = 0.3, alpha = 0.7) + geom_segment(data = fan, aes(x0, y0, xend = x1, yend = y1), colour = "#5d6b61", linewidth = 0.4, alpha = 0.9) + geom_point(data = fan, aes(x1, y1), colour = "#5d6b61", size = 1.3) + geom_segment(data = used_seg, aes(x0, y0, xend = x1, yend = y1), colour = "#b5534e", linewidth = 0.9, arrow = arrow(length = unit(0.16, "cm"), type = "closed")) + geom_point(data = fan, aes(x0, y0), colour = "#16241d", size = 1.8) + coord_equal(xlim = c(0, 100), ylim = c(0, 100), expand = FALSE) + labs(x = "easting", y = "northing") + theme_minimal(base_size = 11) + theme(panel.grid = element_blank(), axis.text = element_text(colour = "#46604a"), legend.title = element_text(size = 8, colour = "#46604a")) ``` ## The conditional likelihood, by hand Within a stratum the animal chose the used step out of the `K + 1` on offer. If selection is proportional to `exp(z * beta)`, the probability that this particular step was the one taken is a softmax over the stratum: $$ \Pr(\text{used}_s) = \frac{\exp(z_{s,\text{used}}\,\beta)}{\sum_{j} \exp(z_{s,j}\,\beta)} . $$ The starting location enters the numerator and every term of the denominator identically, so it cancels. That is the whole point: the local abundance of the resource, which differs from stratum to stratum, drops out, and what remains is the selection for the resource over and above what was locally available. Summing the log of this over strata gives the conditional log-likelihood, which we hand to `optim`. ```{r} #| label: fit clogit_nll <- function(b) { eta <- Z * b -sum(eta[, 1] - log(rowSums(exp(eta)))) } opt <- optim(0, clogit_nll, method = "BFGS", hessian = TRUE, control = list(reltol = 1e-10)) b_hat <- opt$par se_hat <- sqrt(1 / opt$hessian[1, 1]) ci <- b_hat + c(-1, 1) * 1.96 * se_hat rss <- exp(b_hat) ``` The estimate is `r round(b_hat, 3)` with a standard error of `r round(se_hat, 3)` and a 95% interval of [`r round(ci[1], 2)`, `r round(ci[2], 2)`]. Exponentiating gives the relative selection strength: `exp(beta)` is `r round(rss, 2)`, so a step ending in a cell one standard deviation richer in the resource is about `r round(rss, 1)` times as likely as an otherwise identical step ending in an average cell. This estimator is a conditional logistic regression; equivalently it is a Cox proportional hazards model stratified by step, which is how most software fits it. ## Why the strata are not optional It is tempting to throw the used and available steps into one ordinary logistic regression with a single intercept and forget the pairing. That intercept has to stand for the average availability across the whole track, but availability is not constant: as the animal moves through richer and poorer parts of the surface, the resource near its current position rises and falls. A single intercept cannot follow it, and the selection coefficient absorbs the damage. ```{r} #| label: pooled long_z <- as.vector(t(Z)) long_y <- rep(c(1, rep(0, K)), n_str) g_pool <- glm(long_y ~ long_z, family = binomial) b_pool <- unname(coef(g_pool)[2]); se_pool <- summary(g_pool)$coefficients[2, 2] ``` Pooling collapses the coefficient to `r round(b_pool, 3)`, roughly a fifth of the conditional estimate and close to no selection at all. The signal has not weakened; it has been averaged away by an intercept that could not track the moving baseline. Conditioning on the stratum is what protects it. ## Recovering the truth, and where the bias comes from Because this is a simulation we know the answer, so we can check the method directly. If we draw the available steps from the movement kernel that actually generated the track, rather than from a kernel fitted to the observed steps, the estimate should land on the true value. ```{r} #| label: recovery set.seed(424242) Zctrl_true <- matrix(NA_real_, n_str, K) for (s in seq_along(use_idx)) { i <- use_idx[s] Ls <- rgamma(K, shape = sh_move, scale = sc_move) as <- rvm(K, 0, kap_move) ang <- heading[i] + as Zctrl_true[s, ] <- resource(pos[i, 1] + Ls * cos(ang), pos[i, 2] + Ls * sin(ang)) } Zt <- cbind(Zcase, Zctrl_true) clogit_nll_t <- function(b) { eta <- Zt * b; -sum(eta[, 1] - log(rowSums(exp(eta)))) } opt_t <- optim(0, clogit_nll_t, method = "BFGS", hessian = TRUE, control = list(reltol = 1e-10)) b_hat_true <- opt_t$par; se_hat_true <- sqrt(1 / opt_t$hessian[1, 1]) ``` With the true movement kernel supplying availability the estimate is `r round(b_hat_true, 3)`, sitting on the generative value of `r beta_true`. The fitted-kernel estimate above (`r round(b_hat, 2)`) is a shade lower, and the reason is instructive: the tentative gamma and von Mises were fitted to the observed steps, but the observed steps are already the product of selection, so the availability distribution they imply is slightly wrong. Estimating the movement kernel and the habitat selection together, rather than fixing one and inferring the other, removes this bias. That is integrated step selection analysis, and it is the [next post](../integrated-step-selection-issa/). ```{r} #| label: fig-coef #| fig-cap: "Left: the resource at used steps sits higher than at available steps, the selection signal. Right: the conditional estimate recovers the true coefficient (dashed) when availability comes from the movement kernel, is slightly low with a fitted kernel, and collapses toward zero when the strata are pooled into one logistic regression." #| fig-alt: "Two panels. On the left, two overlapping density curves: used steps shifted right of available steps on the resource axis. On the right, three coefficient estimates with error bars against a dashed vertical line at the true value; the two conditional estimates sit near the line while the pooled estimate sits far to the left near zero." #| fig-width: 6.8 #| fig-height: 3.0 dens_df <- rbind(data.frame(kind = "used steps", z = Zcase), data.frame(kind = "available steps", z = as.vector(Zctrl))) pa <- ggplot(dens_df, aes(z, fill = kind, colour = kind)) + geom_density(alpha = 0.45, linewidth = 0.5) + scale_fill_manual(values = c("used steps" = "#275139", "available steps" = "#93a87f")) + scale_colour_manual(values = c("used steps" = "#275139", "available steps" = "#93a87f")) + labs(x = "resource at step endpoint (SD)", y = "density", fill = NULL, colour = NULL) + theme_minimal(base_size = 11) + theme(panel.grid.minor = element_blank(), legend.position = c(0.72, 0.85), legend.background = element_rect(fill = "#ffffff", colour = NA)) coef_df <- data.frame( model = factor(c("conditional: movement kernel", "conditional: fitted kernel", "pooled: strata ignored"), levels = c("pooled: strata ignored", "conditional: fitted kernel", "conditional: movement kernel")), est = c(b_hat_true, b_hat, b_pool), se = c(se_hat_true, se_hat, se_pool)) pb <- ggplot(coef_df, aes(est, model)) + geom_vline(xintercept = beta_true, linetype = "dashed", colour = "#b5534e") + geom_errorbarh(aes(xmin = est - 1.96 * se, xmax = est + 1.96 * se), height = 0.16, colour = "#275139", linewidth = 0.6) + geom_point(colour = "#275139", size = 2.6) + annotate("text", x = beta_true, y = 3.44, label = "true beta", colour = "#b5534e", size = 3, hjust = -0.1) + labs(x = "selection coefficient", y = NULL) + theme_minimal(base_size = 11) + theme(panel.grid.minor = element_blank()) grid.newpage() pushViewport(viewport(layout = grid.layout(1, 2, widths = unit(c(1, 1), "null")))) print(pa, vp = viewport(layout.pos.row = 1, layout.pos.col = 1)) print(pb, vp = viewport(layout.pos.row = 1, layout.pos.col = 2)) ``` ## Selection at which scale A step selection function and a resource selection function answer different questions, and the numbers show it. Fitting a plain resource selection function here, with available points drawn from across the whole domain rather than from reachable steps, gives a coefficient of a different size again. ```{r} #| label: rsf set.seed(999) used_pts <- trk[-1, ] n_av <- nrow(used_pts) * K av_pts <- data.frame(x = runif(n_av, 0, 100), y = runif(n_av, 0, 100)) rsf_y <- c(rep(1, nrow(used_pts)), rep(0, n_av)) rsf_z <- c(resource(used_pts$x, used_pts$y), resource(av_pts$x, av_pts$y)) g_rsf <- glm(rsf_y ~ rsf_z, family = binomial) b_rsf <- unname(coef(g_rsf)[2]) ``` The domain-wide resource selection function returns `r round(b_rsf, 2)`. It is not wrong, and it is not the same thing: it measures broad-scale use of the landscape, the tendency to end up in rich areas over the whole track, whereas the step selection coefficient measures the choice made at each move given where the animal already was. These are different orders of selection in the sense of Johnson (1980), and a study should be clear about which one it is asking for. The step scale is where movement and habitat interact, and it is where the [next post](../integrated-step-selection-issa/) puts them together. ## References Fortin D, Beyer HL, Boyce MS, Smith DW, Duchesne T, Mao JS 2005. Ecology 86(5):1320-1330 (10.1890/04-0953). Thurfjell H, Ciuti S, Boyce MS 2014. Movement Ecology 2:4 (10.1186/s40462-014-0004-7). Forester JD, Im HK, Rathouz PJ 2009. Ecology 90(12):3554-3565 (10.1890/08-0874.1). Johnson DH 1980. Ecology 61(1):65-71 (10.2307/1937156). Manly BFJ, McDonald LL, Thomas DL, McDonald TL, Erickson WP 2002. Resource Selection by Animals, 2nd edn. Kluwer (ISBN 978-1-4020-0677-7). Best DJ, Fisher NI 1979. Journal of the Royal Statistical Society Series C 28(2):152-157 (10.2307/2346732). ## Related tutorials - [Resource selection functions](../resource-selection-functions/) - [Step lengths and turning angles](../step-lengths-turning-angles/) - [Correlated random walks and movement](../correlated-random-walk/) - [Integrated step selection analysis](../integrated-step-selection-issa/)