Practical 4: Multispecies models
Olivier Gimenez (adapted from Marco Salvatori’s code)
May 2023
1 Objective
- Fit occupancy models with co-occurrences between 2 or more species
2 Pre-requisites
Load some packages:
library(tidyverse) # manipulate and visualize data
theme_set(theme_light(base_size = 16)) # set theme for visualisation and font size
library(unmarked) # occupancy modelling3 Data
3.1 Description
We use data from the paper Co-occurrence of snow leopard, wolf and Siberian ibex under livestock encroachment into protected areas across the Mongolian Altai by Salvatori et al. published in Biological Conservation and available at https://www.sciencedirect.com/science/article/abs/pii/S0006320721003463. The authors used camera-traps to sample four areas with different protection status across the Mongolian Altai Mountains, and targeted a predator-prey system composed of livestock, one large herbivore, the Siberian ibex, and two large carnivores, the snow leopard and the wolf. Here we will focus on livestock and wolf. They showed attractino between livestock and wolf, with the potential for human-wildlife conflicts.
We thank Marco Salvatori for sharing the data.
The authors used a two-step approach to infer co-occurrence between species. In the first step they selected the most supported covariates for each species using single-species occupancy models. They first tested the distance to the closest herders’ settlement and the camera performance as covariates on detection probability. They then tested the distance to the closest settlement, terrain slope, elevation, and the sampling area as categorical variable, on occurrence probability while keeping the best encounter structure previously selected. In the second step, the most supported covariates from single-species occupancy models were included as predictors of occurrence and detectability of single species in the multi-species model.
3.2 Read in data
Read in covariates data:
cov <- read.csv2('data/covs_original_ALL.csv', stringsAsFactors = F)
cov <- cov[,c('Sampling.Unit', 'Area', 'CT.sens', 'EL', 'SL', 'Dis')]We have the site id (Sampling.Unit), area, camera performance (CT.sens), elevation (EL), terrain slope (SL) and distance to closest settlement (Dis).
Standardize using decostand() function in
vegan package:
cov[,4:6] <- apply(X = cov[,4:6], MARGIN = 2, FUN = vegan::decostand, method = 'standardize')Read in detections and non-detections:
y <- readRDS(file = 'data/matrici_complessive.rds')Select data for livestock and wolf
y_livestock <- y$domestici
y_wolf <- y$lupo4 Single species analyses
4.1 Organise data
df_livestock <- unmarkedFrameOccu(y = y_livestock, siteCovs = cov)
df_wolf <- unmarkedFrameOccu(y = y_wolf, siteCovs = cov)4.2 Livestock covariates selection
Start with detection:
mod0 <- occu(~1~Area+EL+SL+Dis, df_livestock)
mod1 <- occu(~Dis~Area+EL+SL+Dis, df_livestock)
mod2 <- occu(~CT.sens~Area+EL+SL+Dis, df_livestock)
mod3 <- occu(~Dis+CT.sens~Area+EL+SL+Dis, df_livestock)Gather models and figure out best one:
mods <- fitList('Zero' = mod0,
'Dis' = mod1,
'CT.sens' = mod2,
'Dis + CT.sens' = mod3)
modSel(mods)## Hessian is singular.## nPars AIC delta AICwt cumltvWt
## Dis + CT.sens 10 5230.07 0.00 1.0e+00 1.00
## CT.sens 9 5317.00 86.93 1.3e-19 1.00
## Zero 8 5363.11 133.03 1.3e-29 1.00
## Dis 9 5832.17 602.09 1.8e-131 1.00Proceed with occupancy
mod1 <- occu(~Dis+CT.sens~1, df_livestock)
mod2 <- occu(~Dis+CT.sens~Area, df_livestock)
mod3 <- occu(~Dis+CT.sens~Dis, df_livestock)
mod4 <- occu(~Dis+CT.sens~EL, df_livestock)
mod5 <- occu(~Dis+CT.sens~SL, df_livestock)
mod6 <- occu(~Dis+CT.sens~Area+SL, df_livestock)
mod7 <- occu(~Dis+CT.sens~Area+EL, df_livestock)
mod8 <- occu(~Dis+CT.sens~Area+Dis, df_livestock)
mod9 <- occu(~Dis+CT.sens~EL+Dis, df_livestock)
mod10 <- occu(~Dis+CT.sens~EL+SL, df_livestock)
mod11 <- occu(~Dis+CT.sens~Dis+SL, df_livestock)
mod12 <- occu(~Dis+CT.sens~Area+EL+SL, df_livestock)
mod13 <- occu(~Dis+CT.sens~Area+EL+Dis, df_livestock)
mod14 <- occu(~Dis+CT.sens~Area+SL+Dis, df_livestock)
mod15 <- occu(~Dis+CT.sens~SL+EL+Dis, df_livestock)
mod16 <- occu(~Dis+CT.sens~ Area+EL+SL+Dis, df_livestock)Gather models and figure out best one:
mdls<-fitList('Zero'=mod0,
'K'=mod1,
'Area'=mod2,
'Dis'=mod3,
'EL'=mod4,
'SL'=mod5,
'Area+SL'=mod6,
'Area+EL'=mod7,
'Area+Dis'=mod8,
'EL+Dis'=mod9,
'EL+SL'=mod10,
'Dis+SL'=mod11,
'Area+EL+SL'=mod12,
'Area+EL+Dis'=mod13,
'Area+SL+Dis'=mod14,
'SL+EL+Dis'=mod15,
'Area+EL+SL+Dis'=mod16)
ms <- modSel(mdls)Inspect best model:
summary(mod13)##
## Call:
## occu(formula = ~Dis + CT.sens ~ Area + EL + Dis, data = df_livestock)
##
## Occupancy (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) 0.131 0.277 0.471 0.63789
## AreaSB -0.688 0.416 -1.654 0.09803
## AreaSU 1.008 0.475 2.121 0.03389
## AreaTB -0.811 0.459 -1.765 0.07758
## EL -0.582 0.187 -3.111 0.00187
## Dis -0.617 0.190 -3.241 0.00119
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) -2.436 0.0504 -48.34 0.00e+00
## Dis -0.467 0.0552 -8.47 2.53e-17
## CT.sensmedium 0.487 0.0850 5.73 1.02e-08
##
## AIC: 5228.533
## Number of sites: 216
## optim convergence code: 0
## optim iterations: 87
## Bootstrap iterations: 04.3 Wolf covariates selection
Start with detection:
mod0 <- occu(~1~Area+EL+SL+Dis, df_wolf)
mod1 <- occu(~Dis~Area+EL+SL+Dis, df_wolf)
mod2 <- occu(~CT.sens~Area+EL+SL+Dis, df_wolf)
mod3 <- occu(~Dis+CT.sens~Area+EL+SL+Dis, df_wolf)Gather models and determine best one:
mods <- fitList('Zero'=mod0,
'Dis'=mod1,
'CT.sens'=mod2,
'Dis+CT.sens'=mod3)
modSel(mods)## nPars AIC delta AICwt cumltvWt
## Zero 8 1043.98 0.00 0.523 0.52
## CT.sens 9 1045.83 1.85 0.208 0.73
## Dis 9 1045.98 2.00 0.193 0.92
## Dis+CT.sens 10 1047.82 3.84 0.077 1.00Proceed with occupancy:
## psi
mod1 <- occu(~1~1, df_wolf)
mod2 <- occu(~1~Area, df_wolf)
mod3 <- occu(~1~Dis, df_wolf)
mod4 <- occu(~1~EL, df_wolf)
mod5 <- occu(~1~SL, df_wolf)
mod6 <- occu(~1~Area+SL, df_wolf)
mod7 <- occu(~1~Area+EL, df_wolf)
mod8 <- occu(~1~Area+Dis, df_wolf)
mod9 <- occu(~1~EL+Dis, df_wolf)
mod10 <- occu(~1~EL+SL, df_wolf)
mod11 <- occu(~1~Dis+SL, df_wolf)
mod12 <- occu(~1~Area+EL+SL, df_wolf)
mod13 <- occu(~1~Area+EL+Dis, df_wolf)
mod14 <- occu(~1~Area+SL+Dis, df_wolf)
mod15 <- occu(~1~SL+EL+Dis, df_wolf)
mod16 <- occu(~1~ Area+EL+SL+Dis, df_wolf)Gather models, and determine best one:
mdls <- fitList('Zero'=mod0,
'K'=mod1,
'Area'=mod2,
'Dis'=mod3,
'EL'=mod4,
'SL'=mod5,
'Area+SL'=mod6,
'Area+EL'=mod7,
'Area+Dis'=mod8,
'EL+Dis'=mod9,
'EL+SL'=mod10,
'Dis+SL'=mod11,
'Area+EL+SL'=mod12,
'Area+EL+Dis'=mod13,
'Area+SL+Dis'=mod14,
'SL+EL+Dis'=mod15,
'Area+EL+SL+Dis'=mod16)
modSel(mdls)## nPars AIC delta AICwt cumltvWt
## EL+Dis 4 1039.67 0.00 0.2542 0.25
## Dis 3 1039.87 0.20 0.2303 0.48
## SL+EL+Dis 5 1041.52 1.85 0.1007 0.59
## Dis+SL 4 1041.75 2.08 0.0899 0.67
## Area+EL+Dis 7 1042.21 2.54 0.0712 0.75
## K 2 1043.29 3.62 0.0416 0.79
## Area+Dis 6 1043.47 3.80 0.0379 0.83
## Zero 8 1043.98 4.31 0.0295 0.86
## Area+EL+SL+Dis 8 1043.98 4.31 0.0295 0.88
## Area 5 1044.37 4.70 0.0242 0.91
## EL 3 1044.79 5.12 0.0197 0.93
## Area+EL 6 1045.20 5.53 0.0160 0.94
## Area+SL+Dis 7 1045.23 5.56 0.0158 0.96
## SL 3 1045.24 5.57 0.0157 0.98
## Area+SL 6 1046.16 6.48 0.0099 0.99
## EL+SL 4 1046.73 7.06 0.0074 0.99
## Area+EL+SL 7 1046.98 7.31 0.0066 1.00Inspect best model:
summary(mod9)##
## Call:
## occu(formula = ~1 ~ EL + Dis, data = df_wolf)
##
## Occupancy (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) -0.764 0.221 -3.46 0.000535
## EL -0.299 0.205 -1.46 0.143235
## Dis 0.506 0.198 2.55 0.010718
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## -3.95 0.143 -27.6 8.16e-168
##
## AIC: 1039.671
## Number of sites: 216
## optim convergence code: 0
## optim iterations: 49
## Bootstrap iterations: 04.4 Visualisation
4.4.1 Wolf
Fit best model again:
modW <- occu(~1~EL+Dis, df_wolf)Occupancy vs. distance to settlements:
nd <- data.frame('Area' = factor(rep('KS', 216), levels = c('KS', 'SB', 'SU', 'TB')),
'EL' = 0,
'Dis' = cov$Dis,
'SL' = 0)
wo <- predict(modW, type = 'state', newdata = nd)
wo$Distance <- cov$Dis
ggplot(data = wo, aes(x = Distance, y = Predicted)) +
geom_line() +
geom_ribbon(aes(ymin = lower, ymax = upper ), alpha = .15) +
ylim(0,1) +
labs(y = 'Estimated occupancy', x = 'Distance from settlements', title = 'Wolf') +
scale_x_continuous(breaks = seq(min(wo$Dis), max(wo$Dis), length.out = 3),
labels = c("5000", "10000", "15000"))
Occupancy vs. elevation:
nd <- data.frame('Area' = factor(rep('KS', 216), levels = c('KS', 'SB', 'SU', 'TB')),
'Dis' = 0,
'EL' = cov$EL,
'SL' = 0)
wl <- predict(modW, type='state', newdata = nd)
wl$EL <- cov$EL
ggplot(data = wl, aes(x = EL, y = Predicted)) +
geom_line() +
geom_ribbon(aes(ymin = lower, ymax = upper ), alpha = .15) +
ylim(0,1) +
labs(y = 'Estimated occupancy', x = 'Elevation', title = 'Wolf') +
scale_x_continuous(breaks=seq(min(wl$EL),max(wl$EL), length.out = 5),
labels=c("2000", "2400", "2800","3200", "3500"))
4.4.2 Livestock
Fit best model again:
modLIV <- occu(~Dis+CT.sens~Area+EL+Dis, df_livestock)Occupancy vs. area:
nd <- data.frame('Area' = factor(c('SB', 'TB', 'KS', 'SU'), levels =c('KS', 'SB', 'SU', 'TB') ),
'EL' = rep(0,4),
'SL' = rep(0,4),
'Dis' = rep(0,4))
liv <- predict(modLIV, type = 'state', newdata = nd)
liv$Area<- factor(c('SB', 'TB', 'KS', 'SU'), levels=c('SB', 'TB', 'KS', 'SU'))
ggplot(data = liv, aes(x = Area, y = Predicted)) +
geom_point(aes(x = Area, y = Predicted)) +
geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.5) +
ylim(0,1) +
labs(y = 'Estimated occupancy', title = 'Livestock')
Occupancy vs. elevation:
nd <- data.frame('Area' = factor(rep('KS', 216), levels = c('KS', 'SB', 'SU', 'TB')),
'EL' = cov$EL,
'Dis' = 0,
'SL' = 0)
liv <- predict(modLIV, type = 'state', newdata = nd)
liv$Elevation <- cov$EL
ggplot(data = liv, aes(x = Elevation, y = Predicted)) +
geom_line() +
geom_ribbon(aes(ymin = lower, ymax = upper ), alpha = .15) +
ylim(0,1) +
labs(y = 'Estimated occupancy', x = 'Elevation', title = 'Livestock') +
scale_x_continuous(breaks = c(-2.1,-0.95, 0.3,1.56, 2.5),
labels = c("2000", "2400", "2800", '3200', '3500'))
Occupancy vs. distance to settlements:
nd <- data.frame('Area' = factor(rep('KS', 216), levels = c('KS', 'SB', 'SU', 'TB')),
'EL' = 0,
'Dis' = cov$Dis,
'SL' = 0)
liv2 <- predict(modLIV, type='state', newdata = nd)
liv2$Distance <- cov$Dis
ggplot(data = liv2, aes(x = Distance, y = Predicted)) +
geom_line() +
geom_ribbon(aes(ymin = lower, ymax = upper ), alpha = .15) +
ylim(0,1) +
labs(y = 'Estimated occupancy', x = 'Distance from settlements', title = 'Livestock') +
scale_x_continuous(breaks = seq(min(liv2$Distance), max(liv2$Distance), length.out = 3),
labels = c("5000", "10000", "15000"))
5 Co-occurrence between wolf and livestock
5.1 Prep data
Put together the detection/non-detections for livestock and wolf, along with covariates:
twosp_df <- unmarkedFrameOccuMulti(y = list('livestock' = y_livestock,
'wolf' = y_wolf),
siteCovs = cov)Summarize data:
summary(twosp_df)## unmarkedFrame Object
##
## 216 sites
## 2 species: livestock wolf
## Maximum number of observations per site: 120
## Mean number of observations per site:
## livestock: 65.300925925926 wolf: 65.300925925926
## Sites with at least one detection:
## livestock: 110 wolf: 49
## Tabulation of y observations:
## livestock:
## 0 1 <NA>
## 13314 791 11815
## wolf:
## 0 1 <NA>
## 14015 90 11815
##
## Site-level covariates:
## Sampling.Unit Area CT.sens EL SL
## CT-MON-1-01: 1 KS:63 high :152 Min. :-2.32702 Min. :-1.4793
## CT-MON-1-02: 1 SB:48 medium: 64 1st Qu.:-0.72090 1st Qu.:-0.7736
## CT-MON-1-03: 1 SU:61 Median :-0.05372 Median :-0.1644
## CT-MON-1-04: 1 TB:44 Mean : 0.00000 Mean : 0.0000
## CT-MON-1-05: 1 3rd Qu.: 0.78439 3rd Qu.: 0.4423
## CT-MON-1-06: 1 Max. : 2.59451 Max. : 3.3872
## (Other) :210
## Dis
## Min. :-1.2546
## 1st Qu.:-0.7552
## Median :-0.3236
## Mean : 0.0000
## 3rd Qu.: 0.4638
## Max. : 2.9320
## Visualize design:
plot(twosp_df)
5.2 Model with co-occurrence and constant parameters
We start by fitting a model with constant parameters:
m0 <- occuMulti(detformulas = c('~1','~1'), # detection for livestock, then wolf
stateformulas = c('~1', '~1', '~1'), # marginal occupancy for livestock, wolf, and co-occurrence
data = twosp_df) # data
summary(m0)##
## Call:
## occuMulti(detformulas = c("~1", "~1"), stateformulas = c("~1",
## "~1", "~1"), data = twosp_df, maxOrder = 2L)
##
## Occupancy (logit-scale):
## Estimate SE z P(>|z|)
## [livestock] (Intercept) -0.163 0.187 -0.872 0.383308
## [wolf] (Intercept) -1.063 0.307 -3.463 0.000534
## [livestock:wolf] (Intercept) 0.653 0.385 1.697 0.089662
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## [livestock] (Intercept) -2.17 0.0378 -57.4 0.00e+00
## [wolf] (Intercept) -3.96 0.1451 -27.3 5.04e-164
##
## AIC: 6441.804
## Number of sites: 216
## optim convergence code: 0
## optim iterations: 69
## Bootstrap iterations: 0Get estimated marginal occupancy for wolf:
head(predict(m0, 'state', species = 'wolf'))## Predicted SE lower upper
## 1 0.3296308 0.04726041 0.2418734 0.4268813
## 2 0.3296308 0.04726041 0.2418734 0.4268813
## 3 0.3296308 0.04726041 0.2418734 0.4268813
## 4 0.3296308 0.04726041 0.2418734 0.4268813
## 5 0.3296308 0.04726041 0.2418734 0.4268813
## 6 0.3296308 0.04726041 0.2418734 0.4268813Marginal occupancy for livestock:
head(predict(m0, 'state', species = 'livestock'))## Predicted SE lower upper
## 1 0.5124418 0.03470339 0.4437926 0.5774692
## 2 0.5124418 0.03470339 0.4437926 0.5774692
## 3 0.5124418 0.03470339 0.4437926 0.5774692
## 4 0.5124418 0.03470339 0.4437926 0.5774692
## 5 0.5124418 0.03470339 0.4437926 0.5774692
## 6 0.5124418 0.03470339 0.4437926 0.5774692Now the estimated probability of co-occurrence of wolf and livestock:
head(predict(m0, 'state', species = c('wolf', 'livestock')))## Predicted SE lower upper
## 1 0.2044589 0.04011241 0.1435295 0.2939666
## 2 0.2044589 0.04011241 0.1435295 0.2939666
## 3 0.2044589 0.04011241 0.1435295 0.2939666
## 4 0.2044589 0.04011241 0.1435295 0.2939666
## 5 0.2044589 0.04011241 0.1435295 0.2939666
## 6 0.2044589 0.04011241 0.1435295 0.2939666From these three quantities, we may compute the interaction factor \(\eta\):
psiA <- predict(m0, 'state', species = 'livestock')[1,'Predicted']
psiB <- predict(m0, 'state', species = 'wolf')[1,'Predicted']
psiAB <- predict(m0, 'state', species = c('wolf','livestock'))[1,'Predicted']
eta <- psiAB / (psiA * psiB)
eta## [1] 1.210413There seems to be attraction (convergence).
A way to formally assess the existence of co-occurrence is to compare the model with co-occurrence to the model without co-occurrence. We fit the model without co-occurrence:
m1 <- occuMulti(detformulas = c('~1','~1'), # detection for livestock, then wolf
stateformulas = c('~1', '~1', '~0'), # marginal occupancy for livestock, wolf, and independence
data = twosp_df) # dataNow compare AIC:
tog <- fitList(m0, m1)
modSel(tog)## nPars AIC delta AICwt cumltvWt
## m0 5 6441.80 0.00 0.62 0.62
## m1 4 6442.74 0.94 0.38 1.00The model with co-occurrence has a slightly lower AIC value, suggesting co-occurrence happens between wolf and livetock.
What about wolf occupancy conditional on livestock being present:
head(predict(m0, 'state', species = 'wolf', cond = 'livestock'))## Predicted SE lower upper
## 1 0.3989895 0.0639323 0.2953141 0.5092315
## 2 0.3989895 0.0639323 0.2953141 0.5092315
## 3 0.3989895 0.0639323 0.2953141 0.5092315
## 4 0.3989895 0.0639323 0.2953141 0.5092315
## 5 0.3989895 0.0639323 0.2953141 0.5092315
## 6 0.3989895 0.0639323 0.2953141 0.5092315And wolf occupancy conditional on livestock being absent:
head(predict(m0, 'state', species = 'wolf', cond = '-livestock'))## Predicted SE lower upper
## 1 0.2567323 0.05397693 0.1681009 0.3720474
## 2 0.2567323 0.05397693 0.1681009 0.3720474
## 3 0.2567323 0.05397693 0.1681009 0.3720474
## 4 0.2567323 0.05397693 0.1681009 0.3720474
## 5 0.2567323 0.05397693 0.1681009 0.3720474
## 6 0.2567323 0.05397693 0.1681009 0.3720474We can estimate the realized occupancy of both species and therefore determine sites where both species are co-occurring. By doing so, you can infer hotspots of conflicts.
First get realized occupancy estimates for both species:
re_liv <- ranef(m0, species = "livestock") # estimate posterior distribution of z
z_liv <- bup(re_liv, stat = "mode") # get mode of posterior distribution
re_wolf <- ranef(m0, species = "wolf") # estimate posterior distribution of z
z_wolf <- bup(re_wolf, stat = "mode") # get mode of posterior distributionThen determine sites where both species occur, and proportion
both_species_present <- (z_liv + z_wolf) == 2
cbind(z_liv, z_wolf, both_species_present) ## z_liv z_wolf both_species_present
## [1,] 1 0 0
## [2,] 1 1 1
## [3,] 1 1 1
## [4,] 1 0 0
## [5,] 1 0 0
## [6,] 1 0 0
## [7,] 1 0 0
## [8,] 1 0 0
## [9,] 1 0 0
## [10,] 0 1 0
## [11,] 0 0 0
## [12,] 0 0 0
## [13,] 0 1 0
## [14,] 0 0 0
## [15,] 0 0 0
## [16,] 0 0 0
## [17,] 0 0 0
## [18,] 0 1 0
## [19,] 1 0 0
## [20,] 0 0 0
## [21,] 0 0 0
## [22,] 0 0 0
## [23,] 0 0 0
## [24,] 1 0 0
## [25,] 0 0 0
## [26,] 0 0 0
## [27,] 0 0 0
## [28,] 1 0 0
## [29,] 1 0 0
## [30,] 1 0 0
## [31,] 1 0 0
## [32,] 0 0 0
## [33,] 0 1 0
## [34,] 0 0 0
## [35,] 0 0 0
## [36,] 0 0 0
## [37,] 1 0 0
## [38,] 0 0 0
## [39,] 1 0 0
## [40,] 0 0 0
## [41,] 1 0 0
## [42,] 1 0 0
## [43,] 1 1 1
## [44,] 0 0 0
## [45,] 0 1 0
## [46,] 1 0 0
## [47,] 0 0 0
## [48,] 0 0 0
## [49,] 0 0 0
## [50,] 1 0 0
## [51,] 1 0 0
## [52,] 0 0 0
## [53,] 0 1 0
## [54,] 1 0 0
## [55,] 0 0 0
## [56,] 1 0 0
## [57,] 1 0 0
## [58,] 0 0 0
## [59,] 0 0 0
## [60,] 0 0 0
## [61,] 0 0 0
## [62,] 0 0 0
## [63,] 1 1 1
## [64,] 0 1 0
## [65,] 1 1 1
## [66,] 0 0 0
## [67,] 0 0 0
## [68,] 0 0 0
## [69,] 1 1 1
## [70,] 0 0 0
## [71,] 0 0 0
## [72,] 1 0 0
## [73,] 0 0 0
## [74,] 0 0 0
## [75,] 0 0 0
## [76,] 0 1 0
## [77,] 0 0 0
## [78,] 0 0 0
## [79,] 0 0 0
## [80,] 0 0 0
## [81,] 0 0 0
## [82,] 1 0 0
## [83,] 1 1 1
## [84,] 1 0 0
## [85,] 1 0 0
## [86,] 0 0 0
## [87,] 0 1 0
## [88,] 0 0 0
## [89,] 0 0 0
## [90,] 0 0 0
## [91,] 0 0 0
## [92,] 0 1 0
## [93,] 1 0 0
## [94,] 1 0 0
## [95,] 1 0 0
## [96,] 0 0 0
## [97,] 1 0 0
## [98,] 0 0 0
## [99,] 1 0 0
## [100,] 1 0 0
## [101,] 0 0 0
## [102,] 1 0 0
## [103,] 1 0 0
## [104,] 0 0 0
## [105,] 1 1 1
## [106,] 0 0 0
## [107,] 1 1 1
## [108,] 1 1 1
## [109,] 0 0 0
## [110,] 1 0 0
## [111,] 0 0 0
## [112,] 0 0 0
## [113,] 0 0 0
## [114,] 1 0 0
## [115,] 1 0 0
## [116,] 0 0 0
## [117,] 0 1 0
## [118,] 1 0 0
## [119,] 0 1 0
## [120,] 1 0 0
## [121,] 1 0 0
## [122,] 1 0 0
## [123,] 1 0 0
## [124,] 0 0 0
## [125,] 1 0 0
## [126,] 1 0 0
## [127,] 0 0 0
## [128,] 1 0 0
## [129,] 1 0 0
## [130,] 1 0 0
## [131,] 1 0 0
## [132,] 0 1 0
## [133,] 1 0 0
## [134,] 1 0 0
## [135,] 1 1 1
## [136,] 0 0 0
## [137,] 0 0 0
## [138,] 0 0 0
## [139,] 0 0 0
## [140,] 1 1 1
## [141,] 1 0 0
## [142,] 0 0 0
## [143,] 1 0 0
## [144,] 1 0 0
## [145,] 1 0 0
## [146,] 1 1 1
## [147,] 0 0 0
## [148,] 1 0 0
## [149,] 0 0 0
## [150,] 0 0 0
## [151,] 1 1 1
## [152,] 1 0 0
## [153,] 0 0 0
## [154,] 1 0 0
## [155,] 0 0 0
## [156,] 1 0 0
## [157,] 0 0 0
## [158,] 1 0 0
## [159,] 1 0 0
## [160,] 1 1 1
## [161,] 1 1 1
## [162,] 1 1 1
## [163,] 0 0 0
## [164,] 1 1 1
## [165,] 0 0 0
## [166,] 1 0 0
## [167,] 1 0 0
## [168,] 1 1 1
## [169,] 1 1 1
## [170,] 0 1 0
## [171,] 0 0 0
## [172,] 1 0 0
## [173,] 1 1 1
## [174,] 1 1 1
## [175,] 1 0 0
## [176,] 1 0 0
## [177,] 1 1 1
## [178,] 0 0 0
## [179,] 1 1 1
## [180,] 1 1 1
## [181,] 0 0 0
## [182,] 1 0 0
## [183,] 1 1 1
## [184,] 0 0 0
## [185,] 0 0 0
## [186,] 0 0 0
## [187,] 0 0 0
## [188,] 1 0 0
## [189,] 0 0 0
## [190,] 1 0 0
## [191,] 1 0 0
## [192,] 0 0 0
## [193,] 0 1 0
## [194,] 0 1 0
## [195,] 1 0 0
## [196,] 1 1 1
## [197,] 0 1 0
## [198,] 0 0 0
## [199,] 0 0 0
## [200,] 1 0 0
## [201,] 1 0 0
## [202,] 1 0 0
## [203,] 1 1 1
## [204,] 1 0 0
## [205,] 1 0 0
## [206,] 1 0 0
## [207,] 0 0 0
## [208,] 0 0 0
## [209,] 1 0 0
## [210,] 0 0 0
## [211,] 1 0 0
## [212,] 1 1 1
## [213,] 1 1 1
## [214,] 1 1 1
## [215,] 0 0 0
## [216,] 0 1 0sum(both_species_present) / length(both_species_present) # proportion## [1] 0.14351855.3 Model with co-occurrence and covariates
We fit a model with covariates (see single-species analyses) and co-occurrence between wolf and livestock:
mod_twosp <- occuMulti(detformulas = c('~Dis+CT.sens','~1'),
stateformulas = c('~Area+EL+Dis', '~EL+Dis', '~1'),
data = twosp_df)
summary(mod_twosp)##
## Call:
## occuMulti(detformulas = c("~Dis+CT.sens", "~1"), stateformulas = c("~Area+EL+Dis",
## "~EL+Dis", "~1"), data = twosp_df, maxOrder = 2L)
##
## Occupancy (logit-scale):
## Estimate SE z P(>|z|)
## [livestock] (Intercept) -0.154 0.310 -0.495 0.620390
## [livestock] AreaSB -0.688 0.416 -1.656 0.097648
## [livestock] AreaSU 1.002 0.475 2.111 0.034794
## [livestock] AreaTB -0.813 0.460 -1.768 0.077043
## [livestock] EL -0.548 0.191 -2.873 0.004071
## [livestock] Dis -0.733 0.205 -3.579 0.000345
## [wolf] (Intercept) -1.300 0.341 -3.807 0.000140
## [wolf] EL -0.187 0.209 -0.893 0.372037
## [wolf] Dis 0.599 0.210 2.848 0.004406
## [livestock:wolf] (Intercept) 0.912 0.441 2.068 0.038680
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## [livestock] (Intercept) -2.436 0.0504 -48.35 0.00e+00
## [livestock] Dis -0.467 0.0552 -8.47 2.55e-17
## [livestock] CT.sensmedium 0.487 0.0850 5.72 1.05e-08
## [wolf] (Intercept) -3.939 0.1418 -27.77 1.06e-169
##
## AIC: 6265.704
## Number of sites: 216
## optim convergence code: 0
## optim iterations: 144
## Bootstrap iterations: 0We also fit a model without co-occurrence and compare their respective AIC:
mod_ind <- occuMulti(detformulas = c('~Dis+CT.sens','~1'),
stateformulas = c('~Area+EL+Dis', '~EL+Dis', '~0'), # set interaction to 0
data = twosp_df)
tog <- fitList(mod_twosp, mod_ind)
modSel(tog)## nPars AIC delta AICwt cumltvWt
## mod_twosp 14 6265.70 0.00 0.78 0.78
## mod_ind 13 6268.20 2.50 0.22 1.00The model with co-occurrence has a lower AIC value than the model without co-occurrence.
Now let’s have a look to the estimated wolf occupancy probability conditional on livestock being absent or present:
z01 <- predict(mod_twosp, type = 'state', species = 'wolf', cond = '-livestock', newdata = cov)
z01$speciesB <- 'Conditional on livestock being'
z01$cond <- 'Absent'
z02 <- predict(mod_twosp, type = 'state', species='wolf', cond = 'livestock', newdata = cov)
z02$speciesB <- 'Conditional on livestock being'
z02$cond <- 'Present'
zz2 <- rbind(z01, z02)
zz2 %>% ggplot() +
geom_boxplot(aes(x = cond, y = Predicted)) +
labs(x = "Conditional on livestock being", y = "Estimated wolf occupancy")
6 What if more than 2 species?
I provide data and code for a 3-species example from a paper by Dyck and colleagues entitled Dracula’s ménagerie: A multispecies occupancy analysis of lynx, wildcat, and wolf in the Romanian Carpathians published in Ecology and Evolution available at https://onlinelibrary.wiley.com/doi/full/10.1002/ece3.8921. The code is available at https://github.com/oliviergimenez/occupancy-workshop/tree/main/tutorials/scripts, and the data at https://github.com/oliviergimenez/occupancy-workshop/tree/main/tutorials/data.
