Integrated species distribution models in R

Combine presence-only records with a structured occupancy survey in one likelihood, correcting sampling bias and sharpening the environmental response.

Two kinds of data describe where a species lives, and they have opposite strengths. Opportunistic presence-only records are plentiful but sampled wherever people happen to go, so their density confounds the environment with survey effort. A structured survey visits chosen sites on a fixed protocol, so it is unbiased and carries an absolute scale, but it is expensive and therefore sparse. An integrated model fits both in a single likelihood, letting the structured survey calibrate the bias in the presence-only stream while the presence-only stream lends its coverage to the estimate.

This post builds that model from scratch in base R. We simulate one landscape with a known environmental response, sample it two ways, and compare four estimators: presence-only alone, presence-only with an effort covariate, the survey alone, and the integrated fit.

A landscape and two samples

The intensity of a point process gives the expected number of individuals per unit area. We set a log-linear intensity driven by one environmental covariate, then thin the points by an accessibility surface that stands in for survey effort. Accessibility is correlated with the environment, which is what makes the naive presence-only estimate go wrong.

# a smooth spatial field on a g x g grid, standardised
smooth_field <- function(g, l, seed) {
  set.seed(seed)
  u <- (seq_len(g) - 0.5) / g
  D <- outer(u, u, function(a, b) exp(-(a - b)^2 / (2 * l^2)))
  Sm <- D / rowSums(D)
  as.numeric(scale(as.numeric(Sm %*% matrix(rnorm(g * g), g, g) %*% t(Sm))))
}

g     <- 50                 # grid resolution
Acell <- 1 / (g * g)        # area of one grid cell (unit square)
x     <- smooth_field(g, 0.18, 4711)                 # environmental covariate
raw   <- smooth_field(g, 0.18, 9001)
acc   <- as.numeric(scale(0.6 * x + sqrt(1 - 0.36) * raw))  # accessibility, rho ~ 0.6

b1_true <- 1.0              # true environmental response
a1_true <- 1.6             # effort effect on detection of records
p_true  <- 0.40            # per-visit detection in the survey
Asite   <- 0.0004          # area of one survey site
EN      <- 4000            # expected number of individuals
b0_true <- log(EN / (Acell * sum(exp(b1_true * x))))
round(c(cor_x_acc = cor(x, acc), b0_true = b0_true), 3)

cor_x_acc   b0_true 
    0.630     7.731 

set.seed(1)
lambda <- exp(b0_true + b1_true * x)                 # true intensity surface
N      <- rpois(1, Acell * sum(lambda))
cell   <- sample(length(x), N, replace = TRUE, prob = lambda / sum(lambda))
bias   <- exp(a1_true * acc)                          # effort-driven thinning
keep   <- rbinom(N, 1, pmin(1, bias[cell] / max(bias)))
po     <- tabulate(cell[keep == 1], nbins = length(x))  # presence-only counts per cell

nsite  <- 300
ssite  <- sample(length(x), nsite)                    # structured survey sites
K      <- 4                                            # repeat visits
psi    <- 1 - exp(-Asite * exp(b0_true + b1_true * x[ssite]))  # cloglog occupancy
z      <- rbinom(nsite, 1, psi)
y      <- rbinom(nsite, K, z * p_true)
c(records = sum(po), sites_detected = sum(y > 0), true_occupancy = round(mean(psi), 3))

       records sites_detected true_occupancy 
       288.000        157.000          0.608 

The complementary log-log link is the glue here. If individuals fall as a Poisson process of intensity lambda, the probability that a site of area A holds at least one individual is 1 minus exp(-A times lambda), which is exactly a cloglog occupancy model with an offset of log(A). That shared form lets the survey and the point process speak about the same intensity surface.

Four estimators

We fit each model by maximising its likelihood with optim. The presence-only likelihood is the standard inhomogeneous Poisson process likelihood evaluated on the grid; the occupancy likelihood marginalises over the latent presence state.

# presence-only point process likelihood on the grid
nll_po <- function(par, use_acc) {
  b1 <- par[1]; s0 <- par[2]
  lp <- s0 + b1 * x
  if (use_acc) lp <- lp + par[3] * acc
  mu <- exp(lp)
  -(sum(po * log(mu)) - Acell * sum(mu))
}
fit_po_naive  <- optim(c(1, 0),    nll_po, use_acc = FALSE, method = "BFGS", hessian = TRUE)
fit_po_effort <- optim(c(1, 0, 0), nll_po, use_acc = TRUE,  method = "BFGS", hessian = TRUE)

# survey-only occupancy, marginalising the latent state
nll_occ <- function(par) {
  b0 <- par[1]; b1 <- par[2]; p <- plogis(par[3])
  ps <- 1 - exp(-Asite * exp(b0 + b1 * x[ssite]))
  -sum(log(ps * dbinom(y, K, p) + (1 - ps) * (y == 0)))
}
fit_occ <- optim(c(b0_true, 1, 0), nll_occ, method = "BFGS", hessian = TRUE)

# integrated: shared b0, b1; the records add their own intercept and effort slope
nll_joint <- function(par) {
  b0 <- par[1]; b1 <- par[2]; a0 <- par[3]; p <- plogis(par[4]); a1 <- par[5]
  mu <- exp(b0 + b1 * x + a0 + a1 * acc)
  ll_po  <- sum(po * log(mu)) - Acell * sum(mu)
  ps     <- 1 - exp(-Asite * exp(b0 + b1 * x[ssite]))
  ll_occ <- sum(log(ps * dbinom(y, K, p) + (1 - ps) * (y == 0)))
  -(ll_po + ll_occ)
}
fit_joint <- optim(c(b0_true, 1, -3, 0, 0), nll_joint, method = "BFGS", hessian = TRUE)

se <- function(f, i) sqrt(diag(solve(f$hessian)))[i]
b1_naive  <- fit_po_naive$par[1]
b1_effort <- fit_po_effort$par[1]
b1_occ    <- fit_occ$par[2];   se_occ   <- se(fit_occ, 2)
b1_joint  <- fit_joint$par[2]; se_joint <- se(fit_joint, 2)
round(c(naive = b1_naive, effort = b1_effort, occupancy = b1_occ, integrated = b1_joint), 3)

     naive     effort  occupancy integrated 
     1.389      1.086      1.017      1.078 

The naive presence-only slope is 1.39, well above the truth of 1, because accessibility masquerades as the environment. Add the effort covariate and the point process recovers 1.09, but its intercept mixes the true abundance with the sampling rate, so it has no absolute scale. The survey alone gives 1.02 and does carry a scale, but with a wide standard error of 0.189. The integrated fit matches the survey slope at 1.08 while cutting the standard error to 0.067, about 36 per cent of the survey-only value.

xg <- seq(-2.5, 2.5, length.out = 100)
rel <- function(b1) b1 * xg - b1 * mean(xg)
dl <- rbind(
  data.frame(x = xg, y = rel(b1_true),  who = "truth"),
  data.frame(x = xg, y = rel(b1_naive), who = "naive presence-only"),
  data.frame(x = xg, y = rel(b1_occ),   who = "survey only"),
  data.frame(x = xg, y = rel(b1_joint), who = "integrated"))
ggplot(dl, aes(x, y, colour = who, linetype = who)) +
  geom_line(linewidth = 0.9) +
  scale_colour_manual(values = c("truth" = te$abandoned, "naive presence-only" = te$grazed,
                                 "survey only" = te$sage, "integrated" = te$forest)) +
  scale_linetype_manual(values = c("truth" = "dashed", "naive presence-only" = "solid",
                                   "survey only" = "solid", "integrated" = "solid")) +
  labs(x = "environmental covariate", y = "relative log-intensity",
       title = "The integrated model recovers the true response",
       subtitle = "Naive presence-only overstates the slope; effort correction and integration fix it",
       colour = NULL, linetype = NULL) +
  theme_te() + theme(legend.position = c(0.28, 0.8))

Figure 1: Relative environmental response by estimator. The naive presence-only curve is too steep; the integrated fit tracks the truth.

Does it hold up across data sets?

One realisation could be luck. We repeat the whole exercise many times and record the slope from the survey-only and integrated fits.

run_one <- function(seed) {
  set.seed(seed)
  lam  <- exp(b0_true + b1_true * x)
  Nn   <- rpois(1, Acell * sum(lam))
  cl   <- sample(length(x), Nn, replace = TRUE, prob = lam / sum(lam))
  bi   <- exp(a1_true * acc)
  kp   <- rbinom(Nn, 1, pmin(1, bi[cl] / max(bi)))
  po2  <- tabulate(cl[kp == 1], nbins = length(x))
  ss   <- sample(length(x), nsite)
  ps   <- 1 - exp(-Asite * exp(b0_true + b1_true * x[ss]))
  zz   <- rbinom(nsite, 1, ps); yy <- rbinom(nsite, K, zz * p_true)
  nll_o <- function(par) { b0<-par[1];b1<-par[2];p<-plogis(par[3])
    q<-1-exp(-Asite*exp(b0+b1*x[ss])); -sum(log(q*dbinom(yy,K,p)+(1-q)*(yy==0))) }
  nll_j <- function(par) { b0<-par[1];b1<-par[2];a0<-par[3];p<-plogis(par[4]);a1<-par[5]
    mu<-exp(b0+b1*x+a0+a1*acc); lpo<-sum(po2*log(mu))-Acell*sum(mu)
    q<-1-exp(-Asite*exp(b0+b1*x[ss])); lo<-sum(log(q*dbinom(yy,K,p)+(1-q)*(yy==0))); -(lpo+lo) }
  fo <- tryCatch(optim(c(b0_true,1,0), nll_o, method="BFGS", hessian=TRUE), error=function(e) NULL)
  fj <- tryCatch(optim(c(b0_true,1,-3,0,0), nll_j, method="BFGS", hessian=TRUE), error=function(e) NULL)
  if (is.null(fo) || is.null(fj)) return(rep(NA, 4))
  so <- sqrt(diag(solve(fo$hessian)))[2]; sj <- sqrt(diag(solve(fj$hessian)))[2]
  c(fo$par[2], so, fj$par[2], sj)
}
M  <- 100
mc <- t(sapply(1:M, function(m) run_one(1000 + m)))
mc <- mc[stats::complete.cases(mc), ]
cover <- function(est, s) mean(abs(est - b1_true) <= 1.96 * s)
occ_sd <- sd(mc[, 1]); joint_sd <- sd(mc[, 3])
round(c(occ_mean = mean(mc[, 1]), joint_mean = mean(mc[, 3]),
        occ_sd = occ_sd, joint_sd = joint_sd,
        sd_ratio = joint_sd / occ_sd,
        joint_coverage = cover(mc[, 3], mc[, 4])), 3)

      occ_mean     joint_mean         occ_sd       joint_sd       sd_ratio 
         1.026          1.010          0.171          0.082          0.482 
joint_coverage 
         0.930 

Across 100 data sets both estimators are unbiased, centred near 1, but the integrated slope varies less: its standard deviation is about 48 per cent of the survey-only spread. That is the whole point of integration. You are not choosing between the two data sets; you are letting each cover the other’s weakness.

est <- data.frame(
  method = factor(c("naive\npresence-only", "presence-only\n+ effort", "survey\nonly", "integrated"),
                  levels = c("naive\npresence-only", "presence-only\n+ effort", "survey\nonly", "integrated")),
  b1 = c(b1_naive, b1_effort, b1_occ, b1_joint),
  se = c(se(fit_po_naive, 1), se(fit_po_effort, 1), se_occ, se_joint))
ggplot(est, aes(method, b1)) +
  geom_hline(yintercept = b1_true, colour = te$abandoned, linetype = "dashed") +
  geom_errorbar(aes(ymin = b1 - 1.96 * se, ymax = b1 + 1.96 * se), width = 0.16,
                colour = te$forest, linewidth = 0.7) +
  geom_point(size = 2.8, colour = te$forest) +
  labs(x = NULL, y = "environmental response (b1)",
       title = "Only integration gives both a correct slope and a scale",
       subtitle = "Dashed line is the true response") +
  theme_te()

Figure 2: Slope estimates across four estimators with the truth marked. Presence-only alone is badly biased; the intercept is only identified once the survey enters.

When to reach for this

Integration earns its keep when a large biased data set sits beside a small clean one, and when a plausible covariate explains the bias. It cannot rescue you if the bias is unmodelled and correlated with the environment, and it assumes both streams see the same underlying surface. The next tutorial checks exactly that assumption, showing how a joint model can hide a disagreement between its two data sources.

Related tutorials

• Logistic regression for presence-absence
• GLMs for count data
• Cleaning GBIF occurrence data
• Checking an integrated SDM

References

Dorazio RM (2014) Global Ecology and Biogeography 23: 1472-1484. doi:10.1111/geb.12216

Fithian W, Elith J, Hastie T, Keith DA (2015) Methods in Ecology and Evolution 6: 424-438. doi:10.1111/2041-210X.12242

Koshkina V, Wang Y, Gordon A, Dorazio RM, White M, Stone L (2017) Methods in Ecology and Evolution 8: 420-430. doi:10.1111/2041-210X.12738

Fletcher RJ, Hefley TJ, Robertson EP, Zuckerberg B, McCleery RA, Dorazio RM (2019) Ecology 100: e02710. doi:10.1002/ecy.2710

Miller DAW, Pacifici K, Sanderlin JS, Reich BJ (2019) Methods in Ecology and Evolution 10: 22-37. doi:10.1111/2041-210X.13110

Isaac NJB, Jarzyna MA, Keil P, et al. (2020) Trends in Ecology and Evolution 35: 56-67. doi:10.1016/j.tree.2019.08.006