Skip to contents

Integrated occupancy: combining data sources

Source: vignettes/integrated-occupancy.Rmd
integrated-occupancy.Rmd

Installation 

install.packages("pak")
pak::pak("gcol33/tulpaObs")            # latest from GitHub
pak::pak("gcol33/tulpaObs@v0.1.1")     # a specific tagged release

Tagged releases are listed at https://github.com/gcol33/tulpaObs/releases.

pak resolves the dependency tree, pulling tulpa and tulpaMesh from the r-universe repository DESCRIPTION declares in Additional_repositories, at or above the Imports floor. A C++17 toolchain is needed (Rtools on Windows, Xcode CLI tools on macOS, r-base-dev on Linux); both tulpa and tulpaObs compile their backends on first install.

Two surveys, one species

A structured survey is expensive. A field crew visits a fixed set of sites a fixed number of times, records presence and absence under a known protocol, and the resulting detection histories support a clean occupancy model. The catch is reach: a few hundred sites, often clustered near roads or stations, covering a fraction of the range you care about.

Citizen-science records cover the rest. Birders, naturalists, and app users report sightings across the whole landscape, far beyond where any crew can go. The records are plentiful and broad, but the detection process behind them is murky. Nobody logged the visits that found nothing, effort varies from one observer to the next, and the protocol is whatever each person happened to do that day.

An integrated occupancy model uses both. The two sources disagree about how species are detected, yet they agree about one thing: whether a site is occupied. Occupancy is a property of the place and the species, not of who is looking. So you let the sources share a single latent occupancy state and give each its own detection process. The structured survey anchors the detection calibration; the opportunistic records extend the spatial and covariate reach. int_occu() fits exactly this.

This vignette simulates two sources with a known shared truth, fits the integrated model, reads the per-source detection, and shows where integration sharpens the occupancy estimate over a single source fit alone. It closes with rules for when sources may share a state and when they may not.

The model

Let zi∈{0,1}z_i \in \{0, 1\} be the occupancy state of site ii, shared across every source. It depends on site covariates xix_i through one occupancy regression:

zi∼Bernoulli(ψi),logit(ψi)=xi⊤βocc. z_i \sim \mathrm{Bernoulli}(\psi_i), \qquad \mathrm{logit}(\psi_i) = x_i^\top \beta_{\mathrm{occ}}.

Each source ss then sees that same state through its own detection process. On a visit jj to site ii by source ss, the observation ysijy_{sij} is a detection only when the site is occupied and the source happens to catch it:

ysij∣zi∼Bernoulli(zipsij),logit(psij)=wsij⊤βdet,s. y_{sij} \mid z_i \sim \mathrm{Bernoulli}\!\left(z_i \, p_{sij}\right), \qquad \mathrm{logit}(p_{sij}) = w_{sij}^\top \beta_{\mathrm{det},\,s}.

The occupancy coefficient vector βocc\beta_{\mathrm{occ}} is one object, estimated from all sources at once. The detection coefficients βdet,s\beta_{\mathrm{det},\,s} are a separate vector per source, so a precise survey and a noisy citizen-science stream get their own detection probabilities. A site visited by both sources contributes two detection histories that both inform the one ziz_i. A site seen by only one source contributes through that source alone. The shared ψ\psi is where the pooling happens: every detection history, from every source, is evidence about the same occupancy surface.

This is the structure that makes integration pay. A single sparse source has to estimate ψ\psi and its own pp from a thin set of visits, and the two trade off against each other. Adding a second source that observes overlapping sites brings more visits to bear on zz without forcing the detection processes to match. The occupancy slope sharpens; each psp_s stays free.

The likelihood spells out how the pooling works. For one site, the latent state ziz_i is unobserved, so the model sums over its two values. With every source’s visits to that site collected into the histories y1i,y2i,…y_{1i}, y_{2i}, \dots, the site’s contribution is

ψi∏sPr⁡(ysi∣zi=1)+(1−ψi)𝟏[no detection in any source]. \psi_i \prod_s \Pr(y_{si} \mid z_i = 1) \;+\; (1 - \psi_i)\, \mathbf{1}\!\left[\text{no detection in any source}\right].

The first term is the occupied case, where every source’s detection history enters through its own psp_s. The second is the unoccupied case, possible only when no source detected anything. A single detection in any one source rules out zi=0z_i = 0 for that site, which is why a deeper source can resolve the state of a site that a shallow source left ambiguous. The product over sources inside the occupied term is the channel through which the structured survey’s clean detection sharpens the state estimate that the casual records also depend on.

Simulating two sources

simulate_int_occu() draws a shared occupancy state and then observes it through several sources with their own visit counts and detection probabilities. Here source 1 is a shallow survey with two visits per site; source 2 is a deeper opportunistic stream with five visits. Both observe a block of shared sites, and each also covers sites the other does not.

si <- simulate_int_occu(
  N.total  = 500,                       # unique sites across both sources
  n.data   = 2,                         # two sources
  J        = c(2, 5),                   # visits per source: shallow vs deep
  n.shared = 250,                       # sites observed by both
  beta.occ = c(0.3, 0.9),               # occupancy intercept + slope on x
  beta.det = list(c(0.5, 0), c(0.7, 0)),# per-source detection intercept
  seed     = 11
)

The response si$y is a list with one detection-history matrix per source. They have different visit counts and, in general, different numbers of rows, because each source covers its own set of sites.

str(si$y, max.level = 1)
#> List of 2
#>  $ : int [1:375, 1:2] 1 0 0 1 1 0 0 0 0 1 ...
#>   ..- attr(*, "dimnames")=List of 2
#>  $ : int [1:375, 1:5] 0 0 0 1 1 1 0 0 1 1 ...
#>   ..- attr(*, "dimnames")=List of 2
c(survey_sites = nrow(si$y[[1]]), casual_sites = nrow(si$y[[2]]))
#> survey_sites casual_sites 
#>          375          375

The site covariates live in a single data frame with one row per unique site. si$site_maps records which rows of that frame each source observed, and si$truth carries the occupancy coefficients and the latent state for checking the fit.

str(si$data)
#> 'data.frame':    500 obs. of  1 variable:
#>  $ x: num  -0.591 0.0266 -1.5166 -1.3627 1.1785 ...
si$truth$beta_occ
#> [1] 0.3 0.9

Fitting the integrated model

The call mirrors a single-season occupancy fit. The occupancy formula is the first argument and applies to the shared state. detection is the detection formula. y is the list of source matrices, and naming the list elements names the detection processes in the output. family = int_occu() selects the integrated likelihood.

fit <- tobs(
  ~ x,
  data      = si$data,
  family    = int_occu(),
  detection = ~ 1,
  y         = list(survey = si$y[[1]], casual = si$y[[2]]),
  method    = "laplace",
  control   = list(verbose = FALSE)
)

coef() returns the shared occupancy block plus one detection block per source, each prefixed by the name you gave that source in y.

coef(fit)
#>    psi_(Intercept)              psi_x survey_(Intercept) casual_(Intercept) 
#>          0.1627076          0.7712881          0.4767786          0.7466297

The occupancy slope on x lands near the simulated 0.9. The two detection intercepts come back separately, one labelled survey and one casual, because each source carries its own detection process. summary() adds the posterior standard deviation and credible bounds; the psi_ rows are the shared occupancy predictor, the survey_ and casual_ rows the two detection predictors.

summary(fit)
#>                     estimate  std.error       2.5 %    97.5 %
#> psi_(Intercept)    0.1627076 0.09963058 -0.03256475 0.3579799
#> psi_x              0.7712881 0.10968767  0.55630421 0.9862720
#> survey_(Intercept) 0.4767786 0.10954097  0.26208229 0.6914750
#> casual_(Intercept) 0.7466297 0.06762573  0.61408567 0.8791737

Why integration helps

The payoff is in the occupancy slope. To see it, fit the shallow survey on its own and compare. Source 1 alone is a single-season occupancy model on the sites that source observed, pulled out through site_maps.

rows1 <- si$site_maps[[1]]
d1    <- si$data[rows1, , drop = FALSE]

fit_survey <- tobs(
  ~ x, data = d1, family = occu(),
  detection = ~ 1, y = si$y[[1]],
  method = "laplace", control = list(verbose = FALSE)
)
coef(fit_survey, arm = "psi")
#> (Intercept)           x 
#>   0.3317445   0.7705236

The point estimates are in the same neighbourhood, so the comparison that matters is the uncertainty, and that needs an honest interval. The default Laplace approximation runs its occupancy intervals optimistically narrow on the integrated fit, because the penalised EM weights the pooled data heavily at the mode. For comparing uncertainty across the two fits, run both under NUTS, which samples the full posterior and reports a calibrated spread.

fit_int_nuts <- tobs(
  ~ x, data = si$data, family = int_occu(),
  detection = ~ 1, y = list(survey = si$y[[1]], casual = si$y[[2]]),
  method = "nuts",
  control = list(n.iter = 1000, n.warmup = 500, seed = 1, verbose = FALSE)
)

fit_survey_nuts <- tobs(
  ~ x, data = d1, family = occu(),
  detection = ~ 1, y = si$y[[1]],
  method = "nuts",
  control = list(n.iter = 1000, n.warmup = 500, seed = 1, verbose = FALSE)
)

Read the occupancy slope from each posterior. The integrated fit borrows the deeper source’s visits and the extra sites it covers, so its slope sits closer to the truth of 0.9 and carries a smaller standard deviation and a narrower credible interval.

row_x <- function(f) {
  s <- summary(f)["psi_x", c("estimate", "std.error", "2.5 %", "97.5 %")]
  s <- as.numeric(s)
  c(estimate = s[1], sd = s[2], lower = s[3], upper = s[4],
    width = s[4] - s[3])
}

comparison <- rbind(
  `survey only` = row_x(fit_survey_nuts),
  integrated    = row_x(fit_int_nuts)
)
round(comparison, 3)
#>             estimate    sd lower upper width
#> survey only    0.801 0.169 0.503 1.156 0.653
#> integrated     0.777 0.110 0.574 1.004 0.430

The integrated credible interval is narrower than the survey-only one, and its centre is closer to the true slope. The shallow survey, with two visits per site, struggles to separate occupancy from detection on its own. Folding in the deeper source tightens the shared state without making the two sources pretend to share a detection process.

The gain has a shape worth naming. Integration does not manufacture information; it pools what the sources already carry onto the one quantity they agree about. When the second source is itself thin, or observes few of the same sites, the interval barely moves. The benefit scales with how much real occupancy signal the added source brings to the shared state.

It is worth being precise about what gets sharper and what does not. The occupancy slope tightens because more sites and more visits inform the shared state surface. The detection parameters do not pool across sources, by design, so each source’s detection precision still depends only on that source’s own visits. A casual stream with few visits per site keeps a loose detection estimate even inside a large integrated fit. That separation is the feature, not a shortcoming: the survey’s detection stays calibrated by survey data, and the occupancy estimate borrows strength from both. If you find the casual detection intercept wandering, the fix is more casual visits, not more survey sites.

Source-specific detection

Each source detects the species at its own rate, and the model estimates those rates separately. Make the difference stark: a careful survey with high detection and a casual stream with low detection, observing the same occupancy surface.

sd2 <- simulate_int_occu(
  N.total  = 500, n.data = 2, J = c(5, 4), n.shared = 200,
  beta.occ = c(0.3, 0.9),
  beta.det = list(c(1.2, 0),   # survey: logit p = 1.2, so p about 0.77
                  c(-0.8, 0)), # casual: logit p = -0.8, so p about 0.31
  seed = 21
)

Fit it and read the two detection intercepts back on the probability scale. NUTS keeps the detection blocks honest here.

fit2 <- tobs(
  ~ x, data = sd2$data, family = int_occu(),
  detection = ~ 1, y = list(survey = sd2$y[[1]], casual = sd2$y[[2]]),
  method = "nuts",
  control = list(n.iter = 800, n.warmup = 400, seed = 1, verbose = FALSE)
)

data.frame(
  source       = c("survey", "casual"),
  true_p       = round(plogis(c(1.2, -0.8)), 2),
  fitted_p     = round(plogis(c(coef(fit2, arm = "survey"), coef(fit2, arm = "casual"))), 2),
  row.names    = NULL
)
#>   source true_p fitted_p
#> 1 survey   0.77     0.74
#> 2 casual   0.31     0.34

The casual stream comes back near its true detection of 0.31 and the survey near its higher rate, exactly the contrast the data carry. This is the point of the separate detection blocks: a single pooled detection probability would average the two and misrepresent both. The survey arm can attenuate a little when occupancy and detection trade off along their shared ridge, so read the detection intercepts as a calibration of each source’s reporting rate rather than an exact constant.

The detection formula can also differ in structure across sources, beyond differing in its fitted value. Pass a list of formulas, one per source, to give a source its own detection covariates. A structured survey often records effort or observer, while opportunistic records carry no such fields, so the survey arm can take a covariate the casual arm omits.

sd2$data$effort <- rnorm(nrow(sd2$data))   # a survey-side detection covariate

fit3 <- tobs(
  ~ x, data = sd2$data, family = int_occu(),
  detection = list(survey = ~ effort, casual = ~ 1),
  y = list(survey = sd2$y[[1]], casual = sd2$y[[2]]),
  method = "laplace", control = list(verbose = FALSE)
)
names(coef(fit3))
#> [1] "psi_(Intercept)"    "psi_x"              "survey_(Intercept)"
#> [4] "survey_effort"      "casual_(Intercept)"

The survey block now carries an effort slope alongside its intercept, while the casual block stays intercept-only. Detection covariates must be columns of the site data frame, since the integrated builder reads each source’s detection design from that frame; place any spatial, temporal, or random-effect term on the occupancy formula, where the shared field lives, not on a per-source detection formula.

Practical guidance

Integration is a modelling commitment, not a free lunch. It rests on the assumption that the sources observe the same latent state, and that assumption does real work. A few rules keep it honest.

  • Sources must share the same definition of occupied, over the same closure period. The shared ziz_i is one event: the species is present at site ii during the season both sources sampled. If the structured survey ran in spring and the citizen records span the whole year, they are not measuring the same state, and pooling them biases ψ\psi toward whichever source sampled the wider window. Align the temporal window before you align the model.

  • Overlap is the currency. Aim for a real block of shared sites, not a handful. The deeper source sharpens the shared slope only through sites it observes in common with the others, since those are where two detection histories pin one ziz_i. As a working target, want at least 20 to 30 shared sites, and more when detection is low, so the cross-source calibration has something to stand on. With near-zero overlap the model reduces to fitting each source independently and stacking them, and the integration buys little.

  • Lean on integration when the primary source is thin. The gain is largest when the structured survey is small or shallow, two or three visits per site, a few hundred sites, and a broad opportunistic source can extend its reach. When the structured survey is already deep and wide, a single-source occupancy fit may be cleaner, and the second source mostly adds detection nuisance parameters.

  • Give each source its own detection process, always. The reason to integrate rather than concatenate is that the sources detect differently. Forcing one shared pp across a calibrated survey and a casual stream throws away the calibration the survey provides. Keep the detection blocks separate and let the survey anchor the scale.

  • Read the occupancy slopes, treat the intercepts with care. Occupancy and detection trade off along a ridge, more so with few visits, and the occupancy intercept absorbs part of that trade. The slopes are the inferential payload and recover cleanly; the intercept is a level to report with its interval, not a precise baseline rate. NUTS gives the honest spread when the default Laplace intervals on the integrated fit run narrow.

When not to integrate. The model fails quietly when the sources measure different states. A breeding-season survey and an autumn-migration dataset both record the same species, but presence in one is not presence in the other, and a shared zz forces a contradiction the fit will resolve by splitting the difference. The same applies to a survey that scores a 5×55 \times 5 km cell as occupied and an app that pins a point: the spatial grain differs, so the occupancy units differ, and they cannot share one state without first being brought to a common resolution. When in doubt, ask whether a single coin flip, the one ziz_i, could have generated both sources’ data at that site. If not, fit them separately.

Where to go next