Checking a dispersal kernel

dispersal

seed dispersal

model checking

ecology tutorial

Three checks for a fitted seed dispersal kernel in R: how far the tail rests on a few seeds, whether the trap window truncates it, and if the family matters.

Author

Tidy Ecology

Published

2026-07-23

The three earlier tutorials in this cluster each carried a warning: the ring geometry doubles the scale if you ignore it, the tail family sets long-distance dispersal by fiat, and the mean averages over a tail the data never saw. This one turns those warnings into checks you can run on any fitted kernel. Each targets a specific way a dispersal fit misleads, and each is a few lines of R.

kd_2dt <- function(r, a, p) (2 * p * r / a^2) * (1 + r^2 / a^2)^(-(p + 1))
S_2dt  <- function(R, a, p) (1 + R^2 / a^2)^(-p)
m_2dt  <- function(a, p) if (p > 0.5) p * a * (sqrt(pi) / 2) * gamma(p - 0.5) / gamma(p + 1) else Inf
r2dt   <- function(N, a, p) a * sqrt((1 - runif(N))^(-1 / p) - 1)   # sampler
fit_2dt <- function(r) {
  nll <- function(v) { d <- kd_2dt(r, exp(v[1]), exp(v[2])); -sum(log(pmax(d, 1e-300))) }
  o <- tryCatch(optim(c(log(median(r)), 0), nll, method = "BFGS"),
                error = function(e) optim(c(log(median(r)), 0), nll, method = "Nelder-Mead"))
  list(a = exp(o$par[1]), p = exp(o$par[2])) }

Check 1: how much does the tail rest on a few seeds?

The tail exponent is estimated from the seeds that reached the tail, which are exactly the seeds you have fewest of. A quick way to feel how load-bearing they are: drop the farthest few per cent and refit. If the tail exponent barely moves, it is supported; if it jumps, it was resting on a handful of points.

set.seed(505)
r <- r2dt(500, a = 25, p = 1.3)
f_full <- fit_2dt(r)
keep <- r <= quantile(r, 0.95)          # drop the farthest 5%
f_trim <- fit_2dt(r[keep])
n_drop <- sum(!keep)

Dropping the 25 farthest seeds (five per cent of the sample) moves the tail exponent from 1.33 to 2.83: from a heavy power-law tail to one more than twice as steep. The mean of the fitted kernel falls from 27.5 m to 23.0 m, while the median hardly stirs (20.4 m to 19.5 m). A number that changes this much when a handful of seeds leave the dataset is not something the data have pinned down. The check is not that the seeds are wrong to include; it is that the tail exponent, and everything derived from it, comes with a wide error bar that a single point estimate hides.

Two decreasing curves on a log y-axis against distance. They overlap up to about 60 metres, then the full-data curve stays high while the trimmed curve drops much faster in the tail. — Figure 1: Fitted survival curves before and after dropping the farthest 5% of seeds, on a logarithmic scale. The bulk is unchanged; the tail swings by orders of magnitude on the strength of a few points. The dashed line marks the farthest retained seed.

Check 2: does the trap window cut off the tail?

Seed traps reach only so far. Any seed that lands past the farthest trap is never counted, so the data are truncated at the window edge, and a fit that ignores this treats the truncated sample as if it were the whole kernel. Simulate a seed rain, keep only what falls within a 60 m window, and fit two ways: ignoring the truncation, and with the truncation built into the likelihood.

set.seed(606)
r <- r2dt(800, a = 28, p = 1.2); Dmax <- 60
obs <- r[r <= Dmax]                     # only seeds within the trap window
true_mean <- m_2dt(28, 1.2)
frac_beyond <- mean(r > Dmax); n_beyond <- sum(r > Dmax)

nll_naive <- function(v) -sum(log(pmax(kd_2dt(obs, exp(v[1]), exp(v[2])), 1e-300)))
fn <- optim(c(log(median(obs)), 0), nll_naive, method = "BFGS")
mean_naive <- m_2dt(exp(fn$par[1]), exp(fn$par[2]))

# truncated likelihood: divide the density by the probability of landing inside the window
nll_trunc <- function(v) { a <- exp(v[1]); p <- exp(v[2])
  -sum(log(pmax(kd_2dt(obs, a, p), 1e-300))) + length(obs) * log(1 - S_2dt(Dmax, a, p)) }
ft <- optim(c(log(median(obs)), 0), nll_trunc, method = "BFGS")
a_tr <- exp(ft$par[1]); p_tr <- exp(ft$par[2]); mean_trunc <- m_2dt(a_tr, p_tr)

Of the 800 seeds, 125 (15.6 per cent) landed beyond the 60 m window and never entered the data. The naive fit, treating the truncated sample as complete, reports a mean of 25.1 m against the true 35.1 m, a factor of 1.40 too low. It hits that low number by forcing an implausibly thin tail to fit inside the window. Adding a single truncation term to the likelihood, dividing the density by the probability 1 - S(D_{max}) of landing inside the window, recovers a mean of 38.4 m, back in the right range. The check is to know your window and put it in the model: a fit to trap counts that never mentions the farthest trap distance is underestimating dispersal, and the size of the miss grows with how fat the true tail is.

A histogram of seed distances with the region beyond 60 metres shaded to show it is unobserved. A red curve drops steeply inside the window; a green curve declines more slowly and extends well past 60 metres. — Figure 2: The full seed rain (histogram), with the 60 m trap window shaded. The naive fit clings to the truncated bulk and dies too fast; the truncation-corrected fit recovers a density that carries weight past the window, where 16 per cent of the seeds actually landed.

Check 3: does the answer survive the family?

The final check fits several kernel families to the same data and asks which conclusions hold across all of them. The ones that hold are properties of the data; the ones that swing are choices you made when you picked a family.

kd_gauss <- function(r, a) (2 * r / a^2) * exp(-(r / a)^2)
kd_exp   <- function(r, a) (r / a^2) * exp(-r / a)
kd_ep    <- function(r, a, c) (c * r / (a^2 * gamma(2 / c))) * exp(-(r / a)^c)
S_gauss  <- function(R, a) exp(-(R / a)^2); S_exp <- function(R, a) exp(-R / a) * (1 + R / a)
S_ep     <- function(R, a, c) pgamma((R / a)^c, shape = 2 / c, lower.tail = FALSE)

set.seed(707)
r <- r2dt(550, a = 24, p = 1.6)
fit1 <- function(nll, lo, hi) exp(optimize(nll, c(lo, hi))$minimum)
ag <- fit1(function(la) -sum(log(kd_gauss(r, exp(la)))), -5, 12)
ae <- fit1(function(la) -sum(log(kd_exp(r, exp(la)))), -5, 12)
op <- optim(c(log(median(r)), 0), function(v) -sum(log(pmax(kd_ep(r, exp(v[1]), exp(v[2])), 1e-300))), method = "BFGS")
ap <- exp(op$par[1]); cp <- exp(op$par[2])
ft <- fit_2dt(r); at <- ft$a; pt <- ft$p

# within-data 90th percentile of each fitted kernel, by simulation
set.seed(7071); N <- 2e5
q90 <- c(quantile(ag * sqrt(-log(runif(N))), 0.9), quantile(rgamma(N, 2, scale = ae), 0.9),
         quantile(ap * rgamma(N, 2 / cp, 1)^(1 / cp), 0.9), quantile(r2dt(N, at, pt), 0.9))
# probability of exceeding 150 m under each fitted kernel
p150 <- c(S_gauss(150, ag), S_exp(150, ae), S_ep(150, ap, cp), S_2dt(150, at, pt))
q90_spread <- max(q90) / min(q90)
tail_ratio <- p150[4] / p150[1]

The four families disagree about almost nothing inside the data. Their fitted 90th-percentile distances span 43 m to 46 m, a factor of 1.09, all sitting near the empirical value of 41 m. Move past the data and they part completely: the probability of a seed passing 150 m (beyond the farthest observed seed at 180 m) is 2.7e-03 under the 2Dt and 3.3e-11 under the Gaussian, a factor of 8e+07. A within-data quantile is a stable target that any adequate family agrees on; a far-tail probability, and any long-distance-dispersal claim built on it, is a property of the family you chose, not of the seeds you counted.

A step function and four smooth decreasing curves on a log y-axis against distance. All lines coincide up to about 180 metres; beyond it the four smooth curves fan apart, the 2Dt highest and the Gaussian lowest. — Figure 3: The empirical survival curve (step line) with the four fitted kernels, on a logarithmic scale. Within the data the fits track the steps and each other; past the farthest seed (dashed line) they diverge by orders of magnitude, which is where long-distance-dispersal numbers come from.

Putting the checks together

None of these three checks needs anything beyond the fit you already have. Drop the farthest few per cent and see whether the tail exponent holds. Compare your farthest trap to the spread of the data and, if seeds are piling up near the edge, put a truncation term in the likelihood. Fit two or three families and report only what they agree on, treating the far tail as a modelled scenario rather than a measurement. Run together, they separate the part of a dispersal kernel the data actually support, the bulk and the within-data quantiles, from the part that is extrapolation wearing the clothes of an estimate, the tail and the mean and every long-distance number that leans on them.

References

Nathan R, Klein E, Robledo-Arnuncio JJ, Revilla E 2012. In: Clobert J et al. (eds) Dispersal Ecology and Evolution, Oxford University Press, pp 187-210. ISBN 978-0-19-960889-8.

Clark JS, Silman M, Kern R, Macklin E, HilleRisLambers J 1999. Ecology 80(5):1475-1494.

Bullock JM, Mallada Gonzalez L, Tamme R, Gotzenberger L, White SM, Partel M, Hooftman DAP 2017. Journal of Ecology 105(1):6-19 (10.1111/1365-2745.12666).

--- title: "Checking a dispersal kernel" description: "Three checks for a fitted seed dispersal kernel in R: how far the tail rests on a few seeds, whether the trap window truncates it, and if the family matters." date: "2026-07-23 10:00" categories: [R, dispersal, seed dispersal, model checking, ecology tutorial] image: thumbnail.png image-alt: "An empirical dispersal survival curve with several fitted kernels tracking it within the data and fanning apart in the tail beyond the last seed." --- The three earlier tutorials in this cluster each carried a warning: the [ring geometry](../fitting-dispersal-kernels/) doubles the scale if you ignore it, the [tail family](../fat-tailed-dispersal-kernels/) sets long-distance dispersal by fiat, and the [mean](../the-mean-dispersal-distance/) averages over a tail the data never saw. This one turns those warnings into checks you can run on any fitted kernel. Each targets a specific way a dispersal fit misleads, and each is a few lines of R. ```{r} #| include: false knitr::opts_chunk$set(warning = FALSE, message = FALSE) ``` ```{r} #| label: setup kd_2dt <- function(r, a, p) (2 * p * r / a^2) * (1 + r^2 / a^2)^(-(p + 1)) S_2dt <- function(R, a, p) (1 + R^2 / a^2)^(-p) m_2dt <- function(a, p) if (p > 0.5) p * a * (sqrt(pi) / 2) * gamma(p - 0.5) / gamma(p + 1) else Inf r2dt <- function(N, a, p) a * sqrt((1 - runif(N))^(-1 / p) - 1) # sampler fit_2dt <- function(r) { nll <- function(v) { d <- kd_2dt(r, exp(v[1]), exp(v[2])); -sum(log(pmax(d, 1e-300))) } o <- tryCatch(optim(c(log(median(r)), 0), nll, method = "BFGS"), error = function(e) optim(c(log(median(r)), 0), nll, method = "Nelder-Mead")) list(a = exp(o$par[1]), p = exp(o$par[2])) } ``` ## Check 1: how much does the tail rest on a few seeds? The tail exponent is estimated from the seeds that reached the tail, which are exactly the seeds you have fewest of. A quick way to feel how load-bearing they are: drop the farthest few per cent and refit. If the tail exponent barely moves, it is supported; if it jumps, it was resting on a handful of points. ```{r} #| label: check1 set.seed(505) r <- r2dt(500, a = 25, p = 1.3) f_full <- fit_2dt(r) keep <- r <= quantile(r, 0.95) # drop the farthest 5% f_trim <- fit_2dt(r[keep]) n_drop <- sum(!keep) ``` Dropping the `r n_drop` farthest seeds (five per cent of the sample) moves the tail exponent from `r sprintf("%.2f", f_full$p)` to `r sprintf("%.2f", f_trim$p)`: from a heavy power-law tail to one more than twice as steep. The mean of the fitted kernel falls from `r sprintf("%.1f", m_2dt(f_full$a, f_full$p))` m to `r sprintf("%.1f", m_2dt(f_trim$a, f_trim$p))` m, while the median hardly stirs (`r sprintf("%.1f", median(r))` m to `r sprintf("%.1f", median(r[keep]))` m). A number that changes this much when a handful of seeds leave the dataset is not something the data have pinned down. The check is not that the seeds are wrong to include; it is that the tail exponent, and everything derived from it, comes with a wide error bar that a single point estimate hides. ```{r} #| label: fig-check1 #| echo: false #| fig-cap: "Fitted survival curves before and after dropping the farthest 5% of seeds, on a logarithmic scale. The bulk is unchanged; the tail swings by orders of magnitude on the strength of a few points. The dashed line marks the farthest retained seed." #| fig-alt: "Two decreasing curves on a log y-axis against distance. They overlap up to about 60 metres, then the full-data curve stays high while the trimmed curve drops much faster in the tail." #| fig-width: 7 #| fig-height: 4.2 library(ggplot2) Rg <- seq(1, 260, length.out = 400) d1 <- data.frame(R = rep(Rg, 2), S = c(S_2dt(Rg, f_full$a, f_full$p), S_2dt(Rg, f_trim$a, f_trim$p)), fit = rep(c("all seeds", "farthest 5% dropped"), each = length(Rg))) ggplot(d1, aes(R, S, colour = fit)) + geom_vline(xintercept = max(r[keep]), linetype = "dashed", colour = "#5d6b61") + geom_line(linewidth = 0.9) + scale_y_log10() + scale_colour_manual(values = c("all seeds" = "#275139", "farthest 5% dropped" = "#b5534e")) + labs(x = "distance R (m)", y = "P(distance > R)", colour = NULL) + theme_minimal(base_size = 12) + theme(legend.position = "top", panel.grid.minor = element_blank()) ``` ## Check 2: does the trap window cut off the tail? Seed traps reach only so far. Any seed that lands past the farthest trap is never counted, so the data are truncated at the window edge, and a fit that ignores this treats the truncated sample as if it were the whole kernel. Simulate a seed rain, keep only what falls within a 60 m window, and fit two ways: ignoring the truncation, and with the truncation built into the likelihood. ```{r} #| label: check2 set.seed(606) r <- r2dt(800, a = 28, p = 1.2); Dmax <- 60 obs <- r[r <= Dmax] # only seeds within the trap window true_mean <- m_2dt(28, 1.2) frac_beyond <- mean(r > Dmax); n_beyond <- sum(r > Dmax) nll_naive <- function(v) -sum(log(pmax(kd_2dt(obs, exp(v[1]), exp(v[2])), 1e-300))) fn <- optim(c(log(median(obs)), 0), nll_naive, method = "BFGS") mean_naive <- m_2dt(exp(fn$par[1]), exp(fn$par[2])) # truncated likelihood: divide the density by the probability of landing inside the window nll_trunc <- function(v) { a <- exp(v[1]); p <- exp(v[2]) -sum(log(pmax(kd_2dt(obs, a, p), 1e-300))) + length(obs) * log(1 - S_2dt(Dmax, a, p)) } ft <- optim(c(log(median(obs)), 0), nll_trunc, method = "BFGS") a_tr <- exp(ft$par[1]); p_tr <- exp(ft$par[2]); mean_trunc <- m_2dt(a_tr, p_tr) ``` Of the 800 seeds, `r n_beyond` (`r sprintf("%.1f", 100 * frac_beyond)` per cent) landed beyond the 60 m window and never entered the data. The naive fit, treating the truncated sample as complete, reports a mean of `r sprintf("%.1f", mean_naive)` m against the true `r sprintf("%.1f", true_mean)` m, a factor of `r sprintf("%.2f", true_mean / mean_naive)` too low. It hits that low number by forcing an implausibly thin tail to fit inside the window. Adding a single truncation term to the likelihood, dividing the density by the probability `1 - S(D_{max})` of landing inside the window, recovers a mean of `r sprintf("%.1f", mean_trunc)` m, back in the right range. The check is to know your window and put it in the model: a fit to trap counts that never mentions the farthest trap distance is underestimating dispersal, and the size of the miss grows with how fat the true tail is. ```{r} #| label: fig-check2 #| echo: false #| fig-cap: "The full seed rain (histogram), with the 60 m trap window shaded. The naive fit clings to the truncated bulk and dies too fast; the truncation-corrected fit recovers a density that carries weight past the window, where 16 per cent of the seeds actually landed." #| fig-alt: "A histogram of seed distances with the region beyond 60 metres shaded to show it is unobserved. A red curve drops steeply inside the window; a green curve declines more slowly and extends well past 60 metres." #| fig-width: 7 #| fig-height: 4.2 xr <- seq(0, quantile(r, 0.99), length.out = 500) d2 <- data.frame(r = rep(xr, 2), dens = c(kd_2dt(xr, exp(fn$par[1]), exp(fn$par[2])), kd_2dt(xr, a_tr, p_tr)), fit = rep(c("naive (window ignored)", "truncation corrected"), each = length(xr))) ggplot() + annotate("rect", xmin = Dmax, xmax = max(xr), ymin = 0, ymax = Inf, fill = "#e7e4d8", alpha = 0.6) + geom_histogram(data = data.frame(r = r), aes(r, after_stat(density)), bins = 40, fill = "#dad9ca", colour = "white") + geom_line(data = d2, aes(r, dens, colour = fit), linewidth = 0.9) + geom_vline(xintercept = Dmax, linetype = "dashed", colour = "#5d6b61") + scale_colour_manual(values = c("naive (window ignored)" = "#b5534e", "truncation corrected" = "#275139")) + labs(x = "distance from source (m)", y = "density", colour = NULL) + theme_minimal(base_size = 12) + theme(legend.position = "top", panel.grid.minor = element_blank()) ``` ## Check 3: does the answer survive the family? The final check fits several kernel families to the same data and asks which conclusions hold across all of them. The ones that hold are properties of the data; the ones that swing are choices you made when you picked a family. ```{r} #| label: check3 kd_gauss <- function(r, a) (2 * r / a^2) * exp(-(r / a)^2) kd_exp <- function(r, a) (r / a^2) * exp(-r / a) kd_ep <- function(r, a, c) (c * r / (a^2 * gamma(2 / c))) * exp(-(r / a)^c) S_gauss <- function(R, a) exp(-(R / a)^2); S_exp <- function(R, a) exp(-R / a) * (1 + R / a) S_ep <- function(R, a, c) pgamma((R / a)^c, shape = 2 / c, lower.tail = FALSE) set.seed(707) r <- r2dt(550, a = 24, p = 1.6) fit1 <- function(nll, lo, hi) exp(optimize(nll, c(lo, hi))$minimum) ag <- fit1(function(la) -sum(log(kd_gauss(r, exp(la)))), -5, 12) ae <- fit1(function(la) -sum(log(kd_exp(r, exp(la)))), -5, 12) op <- optim(c(log(median(r)), 0), function(v) -sum(log(pmax(kd_ep(r, exp(v[1]), exp(v[2])), 1e-300))), method = "BFGS") ap <- exp(op$par[1]); cp <- exp(op$par[2]) ft <- fit_2dt(r); at <- ft$a; pt <- ft$p # within-data 90th percentile of each fitted kernel, by simulation set.seed(7071); N <- 2e5 q90 <- c(quantile(ag * sqrt(-log(runif(N))), 0.9), quantile(rgamma(N, 2, scale = ae), 0.9), quantile(ap * rgamma(N, 2 / cp, 1)^(1 / cp), 0.9), quantile(r2dt(N, at, pt), 0.9)) # probability of exceeding 150 m under each fitted kernel p150 <- c(S_gauss(150, ag), S_exp(150, ae), S_ep(150, ap, cp), S_2dt(150, at, pt)) q90_spread <- max(q90) / min(q90) tail_ratio <- p150[4] / p150[1] ``` The four families disagree about almost nothing inside the data. Their fitted 90th-percentile distances span `r sprintf("%.0f", min(q90))` m to `r sprintf("%.0f", max(q90))` m, a factor of `r sprintf("%.2f", q90_spread)`, all sitting near the empirical value of `r sprintf("%.0f", quantile(r, 0.9))` m. Move past the data and they part completely: the probability of a seed passing 150 m (beyond the farthest observed seed at `r sprintf("%.0f", max(r))` m) is `r sprintf("%.1e", p150[4])` under the 2Dt and `r sprintf("%.1e", p150[1])` under the Gaussian, a factor of `r sprintf("%.0e", tail_ratio)`. A within-data quantile is a stable target that any adequate family agrees on; a far-tail probability, and any long-distance-dispersal claim built on it, is a property of the family you chose, not of the seeds you counted. ```{r} #| label: fig-check3 #| echo: false #| fig-cap: "The empirical survival curve (step line) with the four fitted kernels, on a logarithmic scale. Within the data the fits track the steps and each other; past the farthest seed (dashed line) they diverge by orders of magnitude, which is where long-distance-dispersal numbers come from." #| fig-alt: "A step function and four smooth decreasing curves on a log y-axis against distance. All lines coincide up to about 180 metres; beyond it the four smooth curves fan apart, the 2Dt highest and the Gaussian lowest." #| fig-width: 7 #| fig-height: 4.4 rs <- sort(r); emp <- 1 - (seq_along(rs) - 0.5) / length(rs) Rg <- seq(1, 320, length.out = 400) d3 <- data.frame(R = rep(Rg, 4), S = c(S_gauss(Rg, ag), S_exp(Rg, ae), S_ep(Rg, ap, cp), S_2dt(Rg, at, pt)), fam = rep(c("Gaussian", "exponential", "exponential-power", "2Dt"), each = length(Rg))) d3$fam <- factor(d3$fam, levels = c("Gaussian", "exponential", "exponential-power", "2Dt")) ggplot() + geom_step(data = data.frame(R = rs, S = emp), aes(R, S), colour = "#3a3a3a", linewidth = 0.4) + geom_vline(xintercept = max(r), linetype = "dashed", colour = "#5d6b61") + geom_line(data = d3, aes(R, S, colour = fam), linewidth = 0.8) + scale_y_log10(limits = c(1e-8, 1)) + scale_colour_manual(values = c(Gaussian = "#9aa69d", exponential = "#c9a227", `exponential-power` = "#4f8fbf", `2Dt` = "#275139")) + labs(x = "distance R (m)", y = "P(distance > R)", colour = NULL) + theme_minimal(base_size = 12) + theme(legend.position = "top", panel.grid.minor = element_blank()) ``` ## Putting the checks together None of these three checks needs anything beyond the fit you already have. Drop the farthest few per cent and see whether the tail exponent holds. Compare your farthest trap to the spread of the data and, if seeds are piling up near the edge, put a truncation term in the likelihood. Fit two or three families and report only what they agree on, treating the far tail as a modelled scenario rather than a measurement. Run together, they separate the part of a dispersal kernel the data actually support, the bulk and the within-data quantiles, from the part that is extrapolation wearing the clothes of an estimate, the tail and the mean and every long-distance number that leans on them. ## References Nathan R, Klein E, Robledo-Arnuncio JJ, Revilla E 2012. In: Clobert J et al. (eds) Dispersal Ecology and Evolution, Oxford University Press, pp 187-210. ISBN 978-0-19-960889-8. Clark JS, Silman M, Kern R, Macklin E, HilleRisLambers J 1999. Ecology 80(5):1475-1494. Bullock JM, Mallada Gonzalez L, Tamme R, Gotzenberger L, White SM, Partel M, Hooftman DAP 2017. Journal of Ecology 105(1):6-19 (10.1111/1365-2745.12666). ## Related tutorials - [Fitting dispersal kernels in R](../fitting-dispersal-kernels/) - [Fat-tailed dispersal kernels](../fat-tailed-dispersal-kernels/) - [The mean dispersal distance](../the-mean-dispersal-distance/) - [Checking an extreme value model](../checking-an-extreme-value-model/)