Multispecies occupancy models¶

Ken Kellner¶

Outline¶

  1. Background
  2. Model description
  3. Fit model with JAGS
  4. Fit model with unmarked
  5. Marginal and conditional occupancy
  6. Exercise
  7. Model fitting challenges

Background¶

  • Species interactions can drive species distributions
  • Examples:
    • competition
    • predator-prey
    • parasite-host

How can we model interactions?¶

  • Direct observations of interactions (e.g. predation events)
  • Indirect:
    • Over time
    • Over space

Facilited by proliferation of passive detector data (e.g. cameras)

Interactions over time¶

Activity pattern analysis (e.g. Ridout and Linkie 2009, overlap R package)

Doesn't necessarily test if species are typically found in the same locations

Interactions in space¶

Multispecies occupancy models

Don't necessary test if species are active at the same time

Types of multispecies occupancy models¶

  • 2+ species, no interactions explicitly modeled (e.g. community occupancy models; AHM1 Chap 11)
  • Two species, species interaction factor, sometimes has numerical issues (MacKenzie et al. 2004)
  • Two species, asymmetric interactions (Waddle et al. 2010, Richmond et al. 2010)
  • Models above available in PRESENCE, MARK software
  • 2+ species, symmetric interactions (Rota et al. 2016, AHM2 Chap 8) <- focus of this session

Model description¶

  • Suppose you sample for $S = 2$ species at N sites, $i=1...N$
  • Then the latent state $\textbf{Z}_i$ at site $i$ is a vector of length 2, and there are four possible states:
$$ \textbf{Z}_i = \{11\}, \{10\},\{01\}, or \{00\} $$
$$ \textbf{Z}_i \sim \mathrm{MultivariateBernoulli}(\boldsymbol\psi) $$$$ pr(Z_i) = \{11\} \Longrightarrow \psi_{11} $$$$ pr(Z_i) = \{10\} \Longrightarrow \psi_{10} $$$$ pr(Z_i) = \{01\} \Longrightarrow \psi_{01} $$$$ pr(Z_i) = \{00\} \Longrightarrow \psi_{00} $$$$ \sum{\boldsymbol\psi} = 1$$

How do we get the values of $\psi$?¶

Start with the "natural parameters", $f$¶

$f_1$: Term associated with occupancy of species 1 by itself

$f_2$: Term associated with occupancy of species 2 by itself

$f_{12}$: Interaction term between the two species

If $f_{12}$ is positive, presence of one species increases occupancy probability of the other

If $f_{12}$ is negative, presence of one species decreases occupancy probability of the other

Calculating $\psi$: the multinomial logit link¶

Each element of $\mathbf{\psi}$ is a function of the natural parameters $f$.

$$ \begin{align} \psi_{11} & \propto exp(f_{1} + f_{2} + f_{12}) \\ \psi_{10} & \propto exp(f_{1}) \\ \psi_{01} & \propto exp(f_{2}) \\ \psi_{00} & \propto exp(0) = 1 \end{align} $$

Making sure $\mathbf{\psi}$ sums to 1:

$$ \psi_{11} = \frac{exp(f_{1} + f_{2} + f_{12})}{1 + exp(f_{1}) + exp(f_{2}) + exp(f_{1} + f_{2} + f_{12})} $$$$ \psi_{10} = \frac{exp(f_{1})}{1 + exp(f_{1}) + exp(f_{2}) + exp(f_{1} + f_{2} + f_{12})} $$$$ \psi_{01} = \frac{exp(f_{2})}{1 + exp(f_{1}) + exp(f_{2}) + exp(f_{1} + f_{2} + f_{12})} $$$$ \psi_{00} = \frac{1}{1 + exp(f_{1}) + exp(f_{2}) + exp(f_{1} + f_{2} + f_{12})} $$

Modeling detection¶

  • It's the same as single-species occupancy
  • Each species gets its own separate detection model

For example, detection/nondetection of species 1 at site $i$ and occasion $j$ can be modeled as

$ y_{1,ij} \sim \mathrm{Bernoulli}(p_{1} * z_{1,i}) $

where $z_{1,i}$ is the element 1 of $Z_i$ (i.e., latent presence/absence of species 1 at site $i$) and $p_{1}$ is the species 1 detection probability

Multispecies occupancy model in JAGS¶

Look at the data

In [1]:
library(unmarked)
data(MesoCarnivores)
lapply(MesoCarnivores, head) # look at raw data
$bobcat
A matrix: 6 × 3 of type int
X1X2X3
000
000
000
000
000
000
$coyote
A matrix: 6 × 3 of type int
X1X2X3
000
000
000
000
000
000
$redfox
A matrix: 6 × 3 of type int
X1X2X3
000
000
001
000
000
111
$sitecovs
A data.frame: 6 × 6
Dist_5kmHDens_5kmLatitudeLongitudePeople_siteTrail
<dbl><dbl><dbl><dbl><dbl><int>
10.049.3452580.3899441-0.77239580.8571
20.039.4991970.3899250-0.77239200.0020
30.039.6471730.3899111-0.77239540.3871
40.039.5980660.3899166-0.77239720.0030
50.039.6078250.3899179-0.77240040.0000
60.039.7487910.3899058-0.77240460.4431

Format the data for JAGS¶

For now, we will fit a 2-species model (coyote and red fox) with no covariates

In [3]:
jags_data <- list(
    y_coy = MesoCarnivores$coyote,                                  # coyote matrix
    y_fox = MesoCarnivores$redfox,                                  # red fox matrix
    N = nrow(MesoCarnivores$coyote),                                # num sites
    J = ncol(MesoCarnivores$coyote)                                 # num samples
)

Write the model in JAGS language¶

$$ \textbf{Z}_i \sim \mathrm{MultivariateBernoulli}(\boldsymbol\psi) $$

In BUGS code:

Z[i] ~ dcat(Psi[1:4])           # Takes on values 1-4
  1. Both present {11}
  2. Only coyote {10}
  3. Only fox {01}
  4. Neither {00}

Multinomial logit¶

$$ \begin{align} \psi_{11} & \propto exp(f_{1} + f_{2} + f_{12}) \\ \psi_{10} & \propto exp(f_{1}) \\ \psi_{01} & \propto exp(f_{2}) \\ \psi_{00} & \propto 1 \end{align} $$

In BUGS code:

Psi[1] <- exp(f1 + f2 + f12) # {11}
Psi[2] <- exp(f1)            # {10}
Psi[3] <- exp(f2)            # {01}
Psi[4] <- 1                  # {00}

# Priors on natural parameters
f1 ~ dnorm(0, 0.01)
f2 ~ dnorm(0, 0.01)
f12 ~ dnorm(0, 0.01)

JAGS dcat automatically divides Psi by sum(Psi) to make sure Psi sums to 1

Detection model¶

$$ y_{1,ij} \sim \mathrm{Bernoulli}(z_{1,i} * p_{1}) $$$$ y_{2,ij} \sim \mathrm{Bernoulli}(z_{2,i} * p_{2}) $$

In BUGS code:

for (i in 1:N){
  z1[i] <- (Z[i] == 1) + (Z[i] == 2) # is coyote present?
  z2[i] <- (Z[i] == 1) + (Z[i] == 3) # is red fox present?

  for (j in 1:J){
    y_coy[i, j] ~ dbern(z1[i] * p1)
    y_fox[i, j] ~ dbern(z2[i] * p2)
  }
}
p1 ~ dunif(0, 1)
p2 ~ dunif(0, 1)
In [2]:
cat(file="jags_multispecies_null.txt", "
model {

# State process
for (i in 1:N){
  # latent occupancy state
  Z[i] ~ dcat(Psi[1:4])           # Takes on values 1-4
}

# Occupancy probabilities / multinomial logit
Psi[1] <- exp(f1 + f2 + f12) # {11}
Psi[2] <- exp(f1)            # {10}
Psi[3] <- exp(f2)            # {01}
Psi[4] <- 1                  # {00}

# Priors on natural parameters
f1 ~ dnorm(0, 0.01)
f2 ~ dnorm(0, 0.01)
f12 ~ dnorm(0, 0.01)

# Detection process
for (i in 1:N){

  z1[i] <- (Z[i] == 1) + (Z[i] == 2) # is coyote present?
  z2[i] <- (Z[i] == 1) + (Z[i] == 3) # is red fox present?

  for (j in 1:J){
    y_coy[i, j] ~ dbern(z1[i] * p1)
    y_fox[i, j] ~ dbern(z2[i] * p2)
  }
}

# Priors on p
p1 ~ dunif(0, 1)
p2 ~ dunif(0, 1)

}
")
In [4]:
library(jagsUI)

inits <- function(){
  list(Z = rep(1, jags_data$N))
}

params <- c("f1", "f2", "f12", "p1", "p2")

fit <- jags(data = jags_data, inits = inits, parameters.to.save = params,
            model.file = 'jags_multispecies_null.txt', n.chains = 3, n.iter = 1000,
            n.burnin = 500, n.thin = 1, parallel = TRUE)

Processing function input....... 

Done. 
 
Beginning parallel processing using 3 cores. Console output will be suppressed.

Parallel processing completed.

Calculating statistics....... 

Done.
In [5]:
fit
traceplot(fit)

JAGS output for model 'jags_multispecies_null.txt', generated by jagsUI.
Estimates based on 3 chains of 1000 iterations,
adaptation = 100 iterations (sufficient),
burn-in = 500 iterations and thin rate = 1,
yielding 1500 total samples from the joint posterior. 
MCMC ran in parallel for 0.399 minutes at time 2023-07-18 11:52:54.207778.

             mean     sd     2.5%      50%    97.5% overlap0 f  Rhat n.eff
f1         -0.791  0.080   -0.952   -0.790   -0.629    FALSE 1 1.044    51
f2         -2.321  0.149   -2.601   -2.318   -2.035    FALSE 1 1.006   342
f12         1.094  0.212    0.658    1.090    1.495    FALSE 1 1.010   194
p1          0.417  0.018    0.381    0.417    0.452    FALSE 1 1.031    76
p2          0.435  0.029    0.380    0.434    0.492    FALSE 1 1.008   813
deviance 2851.580 63.369 2736.060 2848.895 2983.615    FALSE 1 1.026    97

Successful convergence based on Rhat values (all < 1.1). 
Rhat is the potential scale reduction factor (at convergence, Rhat=1). 
For each parameter, n.eff is a crude measure of effective sample size. 

overlap0 checks if 0 falls in the parameter's 95% credible interval.
f is the proportion of the posterior with the same sign as the mean;
i.e., our confidence that the parameter is positive or negative.

DIC info: (pD = var(deviance)/2) 
pD = 1968.7 and DIC = 4820.273 
DIC is an estimate of expected predictive error (lower is better).

Model with covariates¶

We now need to index Psi and the natural parameters f by site, i.

for (i in 1:N){
  Z[i] ~ dcat(Psi[i, 1:4])           # Takes on values 1-4

  Psi[i,1] <- exp(f1[i] + f2[i] + f12[i]) # {11}
  Psi[i,2] <- exp(f1[i])                  # {10}
  Psi[i,3] <- exp(f2[i])                  # {01}
  Psi[i,4] <- 1                           # {00}

  # ... truncated

}

Linear predictors on f¶

$$ f_{1,i} = \alpha_{f1} + \beta_{f1} * hdens_i $$$$ f_{2,i} = \alpha_{f2} + \beta_{f2} * hdens_i $$$$ f_{12,i} = \alpha_{f12} $$

In BUGS code:

for (i in 1:N){

  # ... truncated

  f1[i] <- alpha_f1 + beta_f1 * hdens[i]
  f2[i] <- alpha_f2 + beta_f2 * hdens[i]
  f12[i] <- alpha_f12
}

Complete model¶

In [6]:
cat(file="jags_multispecies_covs.txt", "
model {

# State process
for (i in 1:N){
  # latent occupancy state
  Z[i] ~ dcat(Psi[i, 1:4])           # Takes on values 1-4
    
  # Occupancy probabilities / multinomial logit
  Psi[i,1] <- exp(f1[i] + f2[i] + f12[i]) # {11}
  Psi[i,2] <- exp(f1[i])                  # {10}
  Psi[i,3] <- exp(f2[i])                  # {01}
  Psi[i,4] <- 1                           # {00}
  
  # Natural parameters
  f1[i] <- alpha_f1 + beta_f1 * hdens[i]
  f2[i] <- alpha_f2 + beta_f2 * hdens[i]
  f12[i] <- alpha_f12
}

# Priors on state parameters
alpha_f1 ~ dunif(-5, 5)
beta_f1 ~ dnorm(0, 0.01)
alpha_f2 ~ dunif(-5, 5)
beta_f2 ~ dnorm(0, 0.01)
alpha_f12 ~ dunif(-5, 5)


# Detection process
for (i in 1:N){

  z1[i] <- (Z[i] == 1) + (Z[i] == 2) # is coyote present?
  z2[i] <- (Z[i] == 1) + (Z[i] == 3) # is red fox present?

  for (j in 1:J){
    y_coy[i, j] ~ dbern(z1[i] * p1)
    y_fox[i, j] ~ dbern(z2[i] * p2)
  }
}

# Priors on p
p1 ~ dunif(0,1)
p2 ~ dunif(0,1)

}
")
In [7]:
jags_data <- list(
    y_coy = MesoCarnivores$coyote,                                  # coyote matrix
    y_fox = MesoCarnivores$redfox,                                  # red fox matrix
    hdens = as.numeric(scale(MesoCarnivores$sitecovs$HDens_5km)),   # covariate
    N = nrow(MesoCarnivores$coyote),                                # num sites
    J = ncol(MesoCarnivores$coyote)                                 # num samples
)
In [8]:
inits <- function(){
  list(Z = rep(1, jags_data$N))
}

params <- c("alpha_f1", "beta_f1", "alpha_f2", "beta_f2",
            "alpha_f12", "p1", "p2")

fit <- jags(data = jags_data, inits = inits, parameters.to.save = params,
            model.file = 'jags_multispecies_covs.txt', n.chains = 3, n.iter = 100,
            n.burnin = 50, n.thin = 1, parallel = TRUE)
fit

Processing function input....... 

Done. 
 
Beginning parallel processing using 3 cores. Console output will be suppressed.

Parallel processing completed.

Calculating statistics....... 

Done. 

JAGS output for model 'jags_multispecies_covs.txt', generated by jagsUI.
Estimates based on 3 chains of 100 iterations,
adaptation = 100 iterations (sufficient),
burn-in = 50 iterations and thin rate = 1,
yielding 150 total samples from the joint posterior. 
MCMC ran in parallel for 0.285 minutes at time 2023-07-18 11:57:12.867507.

              mean     sd     2.5%      50%    97.5% overlap0    f  Rhat n.eff
alpha_f1    -0.776  0.084   -0.942   -0.776   -0.629    FALSE 1.00 1.171    15
beta_f1      0.006  0.063   -0.121    0.005    0.114     TRUE 0.54 1.025   150
alpha_f2    -2.341  0.171   -2.762   -2.333   -2.066    FALSE 1.00 1.095    45
beta_f2      0.549  0.108    0.371    0.524    0.778    FALSE 1.00 1.090    30
alpha_f12    0.989  0.215    0.601    0.977    1.486    FALSE 1.00 1.121   114
p1           0.421  0.021    0.382    0.419    0.460    FALSE 1.00 1.142    19
p2           0.433  0.037    0.362    0.430    0.497    FALSE 1.00 1.026    80
deviance  2849.594 69.478 2715.173 2850.116 2981.134    FALSE 1.00 1.120    27

**WARNING** Rhat values indicate convergence failure. 
Rhat is the potential scale reduction factor (at convergence, Rhat=1). 
For each parameter, n.eff is a crude measure of effective sample size. 

overlap0 checks if 0 falls in the parameter's 95% credible interval.
f is the proportion of the posterior with the same sign as the mean;
i.e., our confidence that the parameter is positive or negative.

DIC info: (pD = var(deviance)/2) 
pD = 2257 and DIC = 5106.6 
DIC is an estimate of expected predictive error (lower is better).
JAGS output for model 'jags_multispecies_covs.txt', generated by jagsUI.
Estimates based on 3 chains of 10000 iterations,
adaptation = 100 iterations (sufficient),
burn-in = 8000 iterations and thin rate = 5,
yielding 1200 total samples from the joint posterior. 
MCMC ran in parallel for 12.987 minutes at time 2023-07-07 08:57:50.122961.

              mean     sd     2.5%      50%    97.5% overlap0     f  Rhat n.eff
alpha_f1    -0.783  0.083   -0.947   -0.782   -0.617    FALSE 1.000 1.002  1200
beta_f1     -0.001  0.074   -0.153   -0.003    0.150     TRUE 0.508 1.003   621
alpha_f2    -2.399  0.158   -2.720   -2.395   -2.105    FALSE 1.000 1.002   711
beta_f2      0.529  0.092    0.365    0.523    0.728    FALSE 1.000 1.000  1200
alpha_f12    1.075  0.230    0.643    1.066    1.534    FALSE 1.000 1.001   996
deviance  2856.690 62.011 2738.322 2855.443 2975.947    FALSE 1.000 1.001  1200

Successful convergence based on Rhat values (all < 1.1). 
Rhat is the potential scale reduction factor (at convergence, Rhat=1). 
For each parameter, n.eff is a crude measure of effective sample size. 

overlap0 checks if 0 falls in the parameter's 95% credible interval.
f is the proportion of the posterior with the same sign as the mean;
i.e., our confidence that the parameter is positive or negative.

DIC info: (pD = var(deviance)/2) 
pD = 1924.6 and DIC = 4781.304 
DIC is an estimate of expected predictive error (lower is better).

Some key conclusions¶

mean     sd     2.5%      50%    97.5% overlap0     f  Rhat n.eff
beta_f1     -0.001  0.074   -0.153   -0.003    0.150     TRUE 0.508 1.003   621

No evidence housing density affects coyote occupancy (remember, species 1 = coyote)

mean     sd     2.5%      50%    97.5% overlap0     f  Rhat n.eff
beta_f2      0.529  0.092    0.365    0.523    0.728    FALSE 1.000 1.000  1200

Positive relationship between housing density and red fox occupancy

mean     sd     2.5%      50%    97.5% overlap0     f  Rhat n.eff
alpha_f12    1.075  0.230    0.643    1.066    1.534    FALSE 1.000 1.001   996

Positive interaction between coyote and red fox occupancy; that is, if coyote is present occupancy of red fox is higher (and vice-versa)

A note on removing interaction terms¶

$$ \begin{align} \psi_{11} & \propto exp(f_{1} + f_{2} + f_{12}) \\ \psi_{10} & \propto exp(f_{1}) \\ \psi_{01} & \propto exp(f_{2}) \\ \psi_{00} & \propto 1 \end{align} $$

What if we force $f_{12} = 0$?

$$ \begin{align} \psi_{11} & \propto exp(f_{1} + f_{2}) \\ \psi_{10} & \propto exp(f_{1}) \\ \psi_{01} & \propto exp(f_{2}) \\ \psi_{00} & \propto 1 \end{align} $$

This is equivalent to fitting two separate occupancy models. I'll show this later.

Thoughts on using JAGS for this model¶

  • It's easy enough for 2 species
  • Important to keep things organized and write yourself comments (e.g. to keep track of which species is which)
  • It gets very complicated very quickly when # species > 2
  • unmarked makes things easier to scale

Multispecies occupancy model in unmarked¶

We will use the function occuMulti

Example data¶

In [9]:
library(unmarked)
In [10]:
data(MesoCarnivores)
lapply(MesoCarnivores, head) # look at raw data
$bobcat
A matrix: 6 × 3 of type int
X1X2X3
000
000
000
000
000
000
$coyote
A matrix: 6 × 3 of type int
X1X2X3
000
000
000
000
000
000
$redfox
A matrix: 6 × 3 of type int
X1X2X3
000
000
001
000
000
111
$sitecovs
A data.frame: 6 × 6
Dist_5kmHDens_5kmLatitudeLongitudePeople_siteTrail
<dbl><dbl><dbl><dbl><dbl><int>
10.049.3452580.3899441-0.77239580.8571
20.039.4991970.3899250-0.77239200.0020
30.039.6471730.3899111-0.77239540.3871
40.039.5980660.3899166-0.77239720.0030
50.039.6078250.3899179-0.77240040.0000
60.039.7487910.3899058-0.77240460.4431

Formatting the data¶

We will create an unmarkedFrameOccuMulti.

The main difference from unmarkedFrameOccu is that you have >1 matrix of observations: we now have one y per species.

Get some more detail by looking at the help file for unmarkedFrameOccuMulti:

In [11]:
?unmarkedFrameOccuMulti

y list¶

The first three elements of MesoCarnivores are the y-matrices. To keep things simple, we'll just use coyote and red fox.

Notice that the resulting list is named: this is important!

In [12]:
ylist <- MesoCarnivores[c(2:3)]
lapply(ylist, head) # look at first few rows
$coyote
A matrix: 6 × 3 of type int
X1X2X3
000
000
000
000
000
000
$redfox
A matrix: 6 × 3 of type int
X1X2X3
000
000
001
000
000
111

Site covariates¶

The last element of the MesoCarnivores list is a data frame of site covariates:

In [13]:
site_covs <- MesoCarnivores$sitecovs
head(site_covs)
A data.frame: 6 × 6
Dist_5kmHDens_5kmLatitudeLongitudePeople_siteTrail
<dbl><dbl><dbl><dbl><dbl><int>
10.049.3452580.3899441-0.77239580.8571
20.039.4991970.3899250-0.77239200.0020
30.039.6471730.3899111-0.77239540.3871
40.039.5980660.3899166-0.77239720.0030
50.039.6078250.3899179-0.77240040.0000
60.039.7487910.3899058-0.77240460.4431

Construct the unmarked frame¶

In [14]:
umf <- unmarkedFrameOccuMulti(y=ylist, siteCovs=site_covs)
head(umf)

Data frame representation of unmarkedFrame object.
Only showing observation matrix for species 1.
   y.1 y.2 y.3 Dist_5km HDens_5km  Latitude  Longitude People_site Trail
1    0   0   0     0.04  9.345258 0.3899441 -0.7723958       0.857     1
2    0   0   0     0.03  9.499197 0.3899250 -0.7723920       0.002     0
3    0   0   0     0.03  9.647173 0.3899111 -0.7723954       0.387     1
4    0   0   0     0.03  9.598066 0.3899166 -0.7723972       0.003     0
5    0   0   0     0.03  9.607825 0.3899179 -0.7724004       0.000     0
6    0   0   0     0.03  9.748791 0.3899058 -0.7724046       0.443     1
7    0   0   0     0.03  9.715359 0.3899097 -0.7724065       0.000     0
8    0   0   0     0.03  9.685091 0.3899135 -0.7724103       0.000     0
9    0   0   0     0.04  9.051558 0.3899955 -0.7724561       0.002     0
10   0   0   0     0.04  9.031836 0.3899971 -0.7724504       0.000     0
In [15]:
summary(umf)

unmarkedFrame Object

1437 sites
2 species: coyote redfox 
Maximum number of observations per site: 3 
Mean number of observations per site:
coyote: 3  redfox: 3  
Sites with at least one detection:
coyote: 401  redfox: 161  
Tabulation of y observations:
coyote:
   0    1 
3685  626 
redfox:
   0    1 
4054  257 

Site-level covariates:
    Dist_5km         HDens_5km           Latitude        Longitude      
 Min.   :0.00000   Min.   :  0.0000   Min.   :0.3454   Min.   :-0.8539  
 1st Qu.:0.00000   1st Qu.:  0.2073   1st Qu.:0.3567   1st Qu.:-0.8009  
 Median :0.00000   Median :  0.7310   Median :0.3753   Median :-0.7940  
 Mean   :0.01419   Mean   :  8.4532   Mean   :0.3726   Mean   :-0.7969  
 3rd Qu.:0.02000   3rd Qu.:  3.0219   3rd Qu.:0.3863   3rd Qu.:-0.7801  
 Max.   :0.13000   Max.   :186.6694   Max.   :0.3967   Max.   :-0.7690  
  People_site          Trail       
 Min.   :0.00000   Min.   :0.0000  
 1st Qu.:0.00000   1st Qu.:0.0000  
 Median :0.00000   Median :0.0000  
 Mean   :0.05807   Mean   :0.3229  
 3rd Qu.:0.00900   3rd Qu.:1.0000  
 Max.   :2.44100   Max.   :1.0000
In [16]:
plot(umf)

Set up the formulas¶

While occu had a single formula for occupancy, occuMulti requires one formula per natural parameter $f$, organized into a stateformulas vector.

Similarly, each species has its own detection probability formula organized into a detformulas vector.

stateformulas¶

It can be hard to keep track of how many natural parameters there are and what each one represents.

It can be helpful to look at the $f$-design matrix, which is generated by unmarkedFrameOccuMulti.

In [17]:
umf@fDesign
A matrix: 4 × 3 of type dbl
f1[coyote]f2[redfox]f3[coyote:redfox]
psi[11]111
psi[10]100
psi[01]010
psi[00]000

The number and order of the formulas in the stateformulas vector should match the column names of this matrix. Therefore we'll need 3 formulas total. For this model we'll set the 1st-order $f$ parameters to be a function of (standardized) housing density. The interaction term will be intercept-only.

Our stateformulas should look like this:

In [18]:
stateformulas <- c("~scale(HDens_5km)","~scale(HDens_5km)","~1")
stateformulas
  1. '~scale(HDens_5km)'
  2. '~scale(HDens_5km)'
  3. '~1'

Notice that the formulas are strings (each wrapped in ""). This is required and is just working around some R limitations.

detformulas¶

This one is easier, there is just one formula per species, so there should be 2 total.

In [19]:
detformulas <- c("~1","~1")
detformulas
  1. '~1'
  2. '~1'

Run occuMulti¶

In [20]:
?occuMulti

We now have all the pieces we need (unmarkedFrameOccuMulti, stateformulas, detformulas) needed to run a model.

In [21]:
mod_hdens <- occuMulti(detformulas=detformulas, stateformulas=stateformulas, data=umf)
mod_hdens

Call:
occuMulti(detformulas = detformulas, stateformulas = stateformulas, 
    data = umf, maxOrder = 2L)

Occupancy:
                             Estimate     SE         z  P(>|z|)
[coyote] (Intercept)        -0.786370 0.0815  -9.64306 5.26e-22
[coyote] scale(HDens_5km)   -0.000668 0.0691  -0.00966 9.92e-01
[redfox] (Intercept)        -2.391812 0.1598 -14.96684 1.21e-50
[redfox] scale(HDens_5km)    0.504456 0.0813   6.20113 5.61e-10
[coyote:redfox] (Intercept)  1.062921 0.2322   4.57859 4.68e-06

Detection:
                     Estimate     SE     z  P(>|z|)
[coyote] (Intercept)   -0.333 0.0764 -4.36 1.31e-05
[redfox] (Intercept)   -0.251 0.1180 -2.13 3.34e-02

AIC: 4800.859
  • Coyote occupancy is higher than red fox
  • No relationship between coyote occupancy and housing density
  • Positive relationship between red fox and housing density
  • Positive interaction between coyote and fox

Occupancy probabilities¶

To get the expected probability $\boldsymbol\psi$ for each occupancy state at each site, use predict. This gives you the probabilities along with standard errors and a 95% CI.

In [22]:
occ_prob <- predict(mod_hdens, type="state")
lapply(occ_prob, head)

Bootstrapping confidence intervals with 100 samples
$Predicted
A matrix: 6 × 4 of type dbl
psi[11]psi[10]psi[01]psi[00]
0.073323330.27259550.055609200.5984720
0.073497270.27249910.055741310.5982624
0.073664800.27240620.055868560.5980605
0.073609160.27243700.055826300.5981275
0.073620220.27243090.055834700.5981142
0.073780030.27234230.055956090.5979216
$SE
A matrix: 6 × 4 of type dbl
psi[11]psi[10]psi[01]psi[00]
0.0092999540.013721850.0082394490.01640244
0.0093183410.013729880.0082507010.01641179
0.0093361110.013738030.0082615900.01642111
0.0093302040.013735270.0082579690.01641798
0.0093313770.013735820.0082586880.01641860
0.0093483680.013743870.0082691090.01642771
$lower
A matrix: 6 × 4 of type dbl
psi[11]psi[10]psi[01]psi[00]
0.058420060.24927980.041962830.5677498
0.058556640.24919840.042057320.5674990
0.058688170.24911990.042148340.5672575
0.058644500.24914600.042118110.5673377
0.058653170.24914080.042124120.5673218
0.058778640.24906590.042210940.5670914
$upper
A matrix: 6 × 4 of type dbl
psi[11]psi[10]psi[01]psi[00]
0.092316030.29570220.070490450.6261836
0.092529150.29554710.070645060.6260361
0.092734390.29539790.070795350.6258938
0.092666240.29544740.070745450.6259411
0.092679780.29543760.070755370.6259317
0.092875560.29529530.070898720.6257959
In [23]:
head(apply(occ_prob$Predicted, 1, sum))
  1. 1
  2. 1
  3. 1
  4. 1
  5. 1
  6. 1

Marginal occupancy¶

It's often more interesting to look at individual species. For example, you might want to know the marginal occupancy of species 2 (red fox) at a site, which is the probability of red fox occupancy regardless of the occupancy states of other species.

This is calculated by summing all $\psi$ states in which red fox is present:

$$ \psi_{marg,i2} = \psi_{11}+\psi_{01} $$

You can do this by specifying the species argument in predict.

In [24]:
redfox_marginal <- predict(mod_hdens, type="state", species="redfox")
head(redfox_marginal)

Bootstrapping confidence intervals with 100 samples
A data.frame: 6 × 4
PredictedSElowerupper
<dbl><dbl><dbl><dbl>
10.12893250.010328710.10908290.1505733
20.12923860.010350210.10935780.1509186
30.12953340.010371120.10962270.1512512
40.12943550.010364160.10953470.1511407
50.12945490.010365540.10955220.1511627
60.12973610.010385610.10980490.1514799

Plotting marginal occupancy¶

Compare across species with a plot:

In [25]:
coy_marginal <- predict(mod_hdens, type='state', species="coyote") # get coyote
marg_plot_dat <- rbind(redfox_marginal[1,], coy_marginal[1,])
marg_plot_dat$Species <- c("Red fox", "Coyote")
marg_plot_dat

Bootstrapping confidence intervals with 100 samples
A data.frame: 2 × 5
PredictedSElowerupperSpecies
<dbl><dbl><dbl><dbl><chr>
0.12893250.010328710.10908290.1505733Red fox
0.34591890.016331980.31077900.3757660Coyote
In [26]:
plot(1:2, marg_plot_dat$Predicted, ylim=c(0.1,0.4), 
     xlim=c(0.5,2.5), pch=19, cex=1.5, xaxt='n', 
     xlab="", ylab="Marginal occupancy and 95% CI")
axis(1, at=1:2, labels=marg_plot_dat$Species)

# CIs
top <- 0.1
for (i in 1:2){
  segments(i, marg_plot_dat$lower[i], i, marg_plot_dat$upper[i])
  segments(i-top, marg_plot_dat$lower[i], i+top)
  segments(i-top, marg_plot_dat$upper[i], i+top)
}

Conditional occupancy¶

Alternatively, you might want to know the probability of occupancy of one species, conditional on the presence of another.

For example, the probability of red fox occupancy, conditional on coyote presence:

$$ \psi_{i, fox, coy} = \frac{\psi_{i{11}}}{\psi_{i{10}} + \psi_{i{11}}}$$

And on coyote absence:

$$ \psi_{i, fox, -coy} = \frac{\psi_{i{01}}}{\psi_{i{01}} + \psi_{i{00}}}$$

With predict, use the species and cond arguments together for this.

In [27]:
redfox_coy <- predict(mod_hdens, type="state", species="redfox", cond="coyote")
head(redfox_coy)

Bootstrapping confidence intervals with 100 samples
A data.frame: 6 × 4
PredictedSElowerupper
<dbl><dbl><dbl><dbl>
10.21196690.024062260.16893180.2612879
20.21242210.024109830.16927230.2618617
30.21286040.024155750.16960010.2624140
40.21271490.024140490.16949130.2622307
50.21274380.024143520.16951290.2622671
60.21316170.024187390.16982550.2627938

What about conditional on coyote absence?

In [28]:
redfox_nocoy <- predict(mod_hdens, type="state", species="redfox", cond="-coyote")
head(redfox_nocoy)

Bootstrapping confidence intervals with 100 samples
A data.frame: 6 × 4
PredictedSElowerupper
<dbl><dbl><dbl><dbl>
10.085018800.011721260.064399610.1066717
20.085230880.011748020.064573540.1068801
30.085435200.011773910.064741140.1070808
40.085367350.011765300.064685480.1070142
50.085380830.011767010.064696540.1070274
60.085575770.011791780.064856470.1072196

Plotting conditional occupancy¶

You can use this output from predict to generate comparison plots.

In [29]:
plot_data <- rbind(redfox_coy[1,], redfox_nocoy[1,])
plot_data$Coyote_status <- c("Present","Absent")
head(plot_data)
A data.frame: 2 × 5
PredictedSElowerupperCoyote_status
<dbl><dbl><dbl><dbl><chr>
10.21196690.024062260.168931790.2612879Present
20.08501880.011721260.064399610.1066717Absent
In [30]:
plot(1:2, plot_data$Predicted, ylim=c(0, 0.3), 
     xlim=c(0.5,2.5), pch=19, cex=1.5, xaxt='n', 
     xlab="Coyote status", ylab="Fox conditional occupancy and 95% CI")
axis(1, at=1:2, labels=plot_data$Coyote_status)

# CIs
top <- 0.1
for (i in 1:2){
  segments(i, plot_data$lower[i], i, plot_data$upper[i])
  segments(i-top, plot_data$lower[i], i+top)
  segments(i-top, plot_data$upper[i], i+top)
}

Plotting covariate effects¶

To plot the effect of housing density on marginal occupancy, we again use predict.

  1. Generate sequence of possible Hdens_5km values for X-axis.
In [31]:
hdens_range <- range(siteCovs(umf)$HDens_5km)
hdens_range
hdens_seq <- seq(hdens_range[1], hdens_range[2], length.out=100)
  1. 0
  2. 186.6694336
  1. predict occupancy at each value of Hdens_5km along our sequence. Supply the sequence as a data.frame to the newdata argument.
In [32]:
nd <- data.frame(HDens_5km = hdens_seq)
occ_hdens_coy <- predict(mod_hdens, type="state", species="coyote", newdata=nd)
occ_hdens_fox <- predict(mod_hdens, type="state", species="redfox", newdata=nd)
head(occ_hdens_coy)
dim(occ_hdens_coy)

Bootstrapping confidence intervals with 100 samples

Bootstrapping confidence intervals with 100 samples
A data.frame: 6 × 4
PredictedSElowerupper
<dbl><dbl><dbl><dbl>
10.34150350.017793570.31152420.3717624
20.34234960.017558410.31297940.3727728
30.34321800.017369810.31445310.3738033
40.34410900.017229330.31594580.3748543
50.34502300.017138210.31745750.3763669
60.34596040.017097310.31860150.3779619
  1. 100
  2. 4

Make the plot¶

In [33]:
plot(hdens_seq, occ_hdens_coy$Predicted, type = 'l', col = 'red',
     xlab = "Housing density", ylab = "Marginal occupancy", ylim = c(0,0.8))
lines(hdens_seq, occ_hdens_fox$Predicted, col = "blue")
legend("bottomright", col=c("red", "blue"), lty=1, legend=c("coyote", "red fox"))

Removing interaction terms¶

In [34]:
stateformulas <- c("~scale(HDens_5km)", "~scale(HDens_5km)")
mod_noint <- occuMulti(detformulas=detformulas, stateformulas=stateformulas, data=umf, maxOrder=1)
mod_noint

Call:
occuMulti(detformulas = detformulas, stateformulas = stateformulas, 
    data = umf, maxOrder = 1)

Occupancy:
                          Estimate     SE      z  P(>|z|)
[coyote] (Intercept)        -0.631 0.0745  -8.47 2.49e-17
[coyote] scale(HDens_5km)    0.108 0.0628   1.72 8.50e-02
[redfox] (Intercept)        -1.921 0.0998 -19.25 1.52e-82
[redfox] scale(HDens_5km)    0.514 0.0811   6.34 2.27e-10

Detection:
                     Estimate     SE     z  P(>|z|)
[coyote] (Intercept)   -0.332 0.0763 -4.36 1.33e-05
[redfox] (Intercept)   -0.261 0.1195 -2.19 2.87e-02

AIC: 4822.196
In [35]:
umf_coy <- unmarkedFrameOccu(y=ylist$coyote, siteCovs=site_covs)
occu(~1~scale(HDens_5km), umf_coy)

Call:
occu(formula = ~1 ~ scale(HDens_5km), data = umf_coy)

Occupancy:
                 Estimate     SE     z  P(>|z|)
(Intercept)        -0.631 0.0745 -8.47 2.51e-17
scale(HDens_5km)    0.108 0.0628  1.72 8.49e-02

Detection:
 Estimate     SE     z  P(>|z|)
   -0.332 0.0763 -4.36 1.32e-05

AIC: 3245.169

Exercise¶

  • Set up a 3-species model with the MesoCarnivores dataset (i.e., add bobcat)
  • Don't include a 3-species interaction term (i.e., force it to be 0)
  • Either write the JAGS model or fit using unmarked
  • I'll post some hints in a few minutes
  • JAGS Math hint:
$$ \begin{aligned} (1) \psi_{111} & \propto exp(f_1 + f_2 + f_3 + f_{12} + f_{13} + f_{23}) \\ (2) \psi_{110} & \propto exp(f_1 + f_2 + f_{12}) \\ (3) \psi_{101} & \propto exp(f_1 + f_3 + f_{13}) \\ (4) \psi_{011} & \propto exp(f_2 + f_3 + f_{23}) \\ (5) \psi_{100} & \propto exp(f_1) \\ (6) \psi_{010} & \propto exp(f_2) \\ (7) \psi_{001} & \propto exp(f_3) \\ (8) \psi_{000} & \propto 1 \end{aligned} $$
  • unmarked hints: set maxOrder argument to occuMulti to 2.
  • You will need 6 formulas in your stateformulas: for f1, f2, f3, f12, f13, f23

JAGS code¶

In [ ]:
cat(file="jags_multispecies_3.txt", "
model {

# State process
for (i in 1:N){

  # latent occupancy state
  Z[i] ~ dcat(Psi[i, 1:8])           # Takes on values 1-8

  # Occupancy probabilities

  Psi[i, 1] <- exp(f1[i] + f2[i] + f3[i] + f12[i] + f13[i] + f23[i]) # {111}
  Psi[i, 2] <- exp(f1[i] + f2[i] + f12[i])                   # {110}
  Psi[i, 3] <- exp(f1[i] + f3[i] + f13[i])                   # {101}
  Psi[i, 4] <- exp(f2[i] + f3[i] + f23[i])                   # {011}
  Psi[i, 5] <- exp(f1[i])                                    # {100}
  Psi[i, 6] <- exp(f2[i])                                    # {010}
  Psi[i, 7] <- exp(f3[i])                                    # {001}
  Psi[i, 8] <- 1                                             # {000}

  # Natural parameters
  f1[i] <- alpha_f1 + beta_f1 * hdens[i]
  f2[i] <- alpha_f2 + beta_f2 * hdens[i]
  f3[i] <- alpha_f3 + beta_f3 * hdens[i]
  f12[i] <- alpha_f12
  f13[i] <- alpha_f13
  f23[i] <- alpha_f23
}

# State priors
alpha_f1 ~ dunif(-5, 5)
beta_f1 ~ dnorm(0, 0.01)
alpha_f2 ~ dunif(-5, 5)
beta_f2 ~ dnorm(0, 0.01)
alpha_f3 ~ dunif(-5, 5)
beta_f3 ~ dnorm(0, 0.01)
alpha_f12 ~ dunif(-5,5)
alpha_f13 ~ dunif(-5,5)
alpha_f23 ~ dunif(-5,5)

# Detection process
for (i in 1:N){

  z1[i] <- (Z[i] == 1) + (Z[i] == 2) + (Z[i] == 3) + (Z[i] == 5) # is coyote present?
  z2[i] <- (Z[i] == 1) + (Z[i] == 2) + (Z[i] == 4) + (Z[i] == 6) # is red fox present?
  z3[i] <- (Z[i] == 1) + (Z[i] == 3) + (Z[i] == 4) + (Z[i] == 7) # is bobcat present?

  for (j in 1:J){
    y_coy[i, j] ~ dbern(z1[i] * p1)
    y_fox[i, j] ~ dbern(z2[i] * p2)
    y_bob[i, j] ~ dbern(z3[i] * p3)
  }

}

# Detection priors
p1 ~ dunif(0,1)
p2 ~ dunif(0,1)
p3 ~ dunif(0,1)


}
")

unmarked code¶

In [36]:
ylist <- MesoCarnivores[c(1:3)]
umf <- unmarkedFrameOccuMulti(y=ylist, siteCovs=site_covs)
stateformulas <- c("~scale(HDens_5km)","~scale(HDens_5km)","~scale(HDens_5km)","~1","~1","~1")
detformulas <- rep("~1", 3)
(mod_3species <- occuMulti(detformulas=detformulas, stateformulas=stateformulas, data=umf, maxOrder = 2))

Call:
occuMulti(detformulas = detformulas, stateformulas = stateformulas, 
    data = umf, maxOrder = 2)

Occupancy:
                            Estimate     SE      z  P(>|z|)
[bobcat] (Intercept)         -2.5364 0.2842  -8.92 4.53e-19
[bobcat] scale(HDens_5km)    -3.6008 0.9191  -3.92 8.94e-05
[coyote] (Intercept)         -1.3439 0.1470  -9.14 6.04e-20
[coyote] scale(HDens_5km)     0.0884 0.0742   1.19 2.34e-01
[redfox] (Intercept)         -2.2555 0.1566 -14.41 4.65e-47
[redfox] scale(HDens_5km)     0.4392 0.0831   5.29 1.24e-07
[bobcat:coyote] (Intercept)   1.8173 0.2714   6.70 2.16e-11
[bobcat:redfox] (Intercept)  -1.3385 0.3869  -3.46 5.42e-04
[coyote:redfox] (Intercept)   1.3740 0.2597   5.29 1.22e-07

Detection:
                     Estimate     SE     z  P(>|z|)
[bobcat] (Intercept)   -1.139 0.1421 -8.02 1.07e-15
[coyote] (Intercept)   -0.328 0.0759 -4.32 1.53e-05
[redfox] (Intercept)   -0.252 0.1179 -2.14 3.26e-02

AIC: 6530.476

Model fitting challenges¶

  • Multispecies occupancy is fairly complicated
  • You might get poor estimates under certain conditions:
    • Sparse data (many 0s)
    • Boundary estimates (occupancy close to 0 or 1)
    • Few observations where multiple species are detected
    • Separation (perfect correlation with covariate)
  • What do I mean by poor estimates? Very large absolute values and SEs
In [37]:
state_complex <- c(rep("~scale(Dist_5km)+scale(HDens_5km)", 6), 0)
det_complex <- rep("~Trail",3)

mod_complex <- occuMulti(stateformulas=state_complex, detformulas=det_complex, umf)
mod_complex

Call:
occuMulti(detformulas = det_complex, stateformulas = state_complex, 
    data = umf, maxOrder = 3L)

Occupancy:
                                 Estimate     SE      z  P(>|z|)
[bobcat] (Intercept)             -23.0171  5.787 -3.977 6.97e-05
[bobcat] scale(Dist_5km)          -2.4249  0.689 -3.518 4.34e-04
[bobcat] scale(HDens_5km)        -82.3836 19.799 -4.161 3.17e-05
[coyote] (Intercept)              -0.6789  0.225 -3.017 2.55e-03
[coyote] scale(Dist_5km)          -0.0176  0.139 -0.127 8.99e-01
[coyote] scale(HDens_5km)         -0.5534  0.748 -0.740 4.59e-01
[redfox] (Intercept)              -1.3946  0.257 -5.425 5.79e-08
[redfox] scale(Dist_5km)          -0.5293  0.250 -2.114 3.45e-02
[redfox] scale(HDens_5km)          0.2108  0.261  0.808 4.19e-01
[bobcat:coyote] (Intercept)        6.7598  6.390  1.058 2.90e-01
[bobcat:coyote] scale(Dist_5km)    1.6979  0.695  2.444 1.45e-02
[bobcat:coyote] scale(HDens_5km)  17.9202 21.460  0.835 4.04e-01
[bobcat:redfox] (Intercept)       15.3983  3.461  4.449 8.64e-06
[bobcat:redfox] scale(Dist_5km)    0.8836  0.439  2.014 4.40e-02
[bobcat:redfox] scale(HDens_5km)  64.1330 12.375  5.182 2.19e-07
[coyote:redfox] (Intercept)        1.1084  0.363  3.050 2.29e-03
[coyote:redfox] scale(Dist_5km)    0.1149  0.340  0.338 7.35e-01
[coyote:redfox] scale(HDens_5km)   0.9046  0.781  1.159 2.47e-01

Detection:
                     Estimate     SE      z  P(>|z|)
[bobcat] (Intercept)    -2.83 0.1419 -19.91 3.12e-88
[bobcat] Trail           1.74 0.1542  11.26 2.10e-29
[coyote] (Intercept)    -1.96 0.0984 -19.88 5.90e-88
[coyote] Trail           2.17 0.1220  17.75 1.63e-70
[redfox] (Intercept)    -1.59 0.1601  -9.93 2.97e-23
[redfox] Trail           1.78 0.1997   8.93 4.12e-19

AIC: 5958.196

Potential Solutions¶

  • Fewer covariates
  • Fewer species
  • Adjust observation period length if possible
  • Set higher-order interaction terms to 0
  • Penalized likelihood with optimizePenalty function (Clipp et al. 2021)

Penalized likelihood¶

"Bayes-inspired" penalty (Hutchinson et al. 2015):

$$ -\lambda\frac{1}{2}\sum_i{}\theta_i^2 $$

Where $\theta$ is the vector of parameter estimates.

Size of penalty increases as values of $\theta$ increase.

Tradeoff:

  • Introduce small bias in parameter estimates
  • Potentially large reduction in variance
mod_penalty <- occuMulti(stateformulas=state_complex, detformulas=det_complex, umf, penalty=1)
mod_penalty

Bootstraping covariance matrix


Call:
occuMulti(detformulas = det_complex, stateformulas = state_complex, 
    data = umf, penalty = 1, maxOrder = 3L)

Occupancy:
                                 Estimate    SE      z  P(>|z|)
[bobcat] (Intercept)              -1.7810 0.269 -6.627 3.43e-11
[bobcat] scale(Dist_5km)          -1.3143 0.293 -4.491 7.08e-06
[bobcat] scale(HDens_5km)         -2.8200 0.539 -5.230 1.70e-07
[coyote] (Intercept)              -0.6049 0.214 -2.830 4.66e-03
[coyote] scale(Dist_5km)           0.0285 0.106  0.270 7.87e-01
[coyote] scale(HDens_5km)         -1.0908 0.378 -2.887 3.88e-03
[redfox] (Intercept)              -1.5659 0.305 -5.140 2.75e-07
[redfox] scale(Dist_5km)          -0.3068 0.157 -1.958 5.02e-02
[redfox] scale(HDens_5km)          0.4730 0.546  0.867 3.86e-01
[bobcat:coyote] (Intercept)        1.1871 0.434  2.736 6.22e-03
[bobcat:coyote] scale(Dist_5km)    0.9347 0.346  2.705 6.84e-03
[bobcat:coyote] scale(HDens_5km)  -0.3218 1.008 -0.319 7.50e-01
[bobcat:redfox] (Intercept)       -0.8831 0.347 -2.545 1.09e-02
[bobcat:redfox] scale(Dist_5km)    0.0364 0.305  0.119 9.05e-01
[bobcat:redfox] scale(HDens_5km)   2.5609 1.074  2.384 1.71e-02
[coyote:redfox] (Intercept)        1.0001 0.310  3.227 1.25e-03
[coyote:redfox] scale(Dist_5km)    0.0236 0.167  0.141 8.87e-01
[coyote:redfox] scale(HDens_5km)   1.3920 0.345  4.030 5.57e-05

Detection:
                     Estimate     SE      z  P(>|z|)
[bobcat] (Intercept)    -2.44 0.1755 -13.89 7.63e-44
[bobcat] Trail           1.74 0.1966   8.87 7.25e-19
[coyote] (Intercept)    -1.89 0.0926 -20.45 5.64e-93
[coyote] Trail           2.10 0.1179  17.84 3.31e-71
[redfox] (Intercept)    -1.49 0.1782  -8.35 6.64e-17
[redfox] Trail           1.72 0.1714  10.04 9.73e-24

AIC: 6135.555

How to choose penalty value? optimizePenalty() function.

mod_penalty <- optimizePenalty(mod_complex, penalties=c(0.5,1))

Literature Cited¶

Clipp, H.L., Evans, A.L., Kessinger, B.E., Kellner, K. and Rota, C.T. (2021). A penalized likelihood for multispecies occupancy models improves predictions of species interactions. Ecology, p.e03520.

MacKenzie, D. I., Bailey, L. L., & Nichols, J. D. (2004). Investigating species co‐occurrence patterns when species are detected imperfectly. Journal of Animal Ecology, 73(3), 546-555.

Richmond, O.M., Hines, J.E. and Beissinger, S.R. (2010). Two‐species occupancy models: a new parameterization applied to co‐occurrence of secretive rails. Ecological Applications, 20(7), pp.2036-2046.

Ridout, M. S., & Linkie, M. (2009). Estimating overlap of daily activity patterns from camera trap data. Journal of Agricultural, Biological, and Environmental Statistics, 14, 322-337.

Rota, C.T., Ferreira, M.A., Kays, R.W., Forrester, T.D., Kalies, E.L., McShea, W.J., Parsons, A.W. and Millspaugh, J.J. (2016). A multispecies occupancy model for two or more interacting species. Methods in Ecology and Evolution, 7(10), pp.1164-1173.

Waddle, J.H., Dorazio, R.M., Walls, S.C., Rice, K.G., Beauchamp, J., Schuman, M.J. and Mazzotti, F.J. (2010). A new parameterization for estimating co‐occurrence of interacting species. Ecological Applications, 20(5), pp.1467-1475.