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 = "inside", legend.position.inside = 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. Two cautions on reading these figures. The 36 per cent quoted earlier came from one fit and one pair of standard errors, so it is a draw from this same distribution rather than the rule. And the integrated fit’s model-based standard error sits a little below its actual sampling spread, which is why the joint interval covers 0.93 rather than 0.95.

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 with 95 per cent intervals from four estimators, against a dashed line at the truth. The naive presence-only interval sits entirely above the truth; the survey-only interval covers it but is much the widest, and the integrated fit keeps the survey’s position with a far narrower interval.

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.

References

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

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

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

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

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

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

Related tutorials

• Logistic regression for presence-absence
• MaxEnt as a Poisson point process
• Cleaning GBIF occurrence data
• Checking an integrated SDM