A multi-species (species interactions) occupancy model

Marginal vs conditional occupancy

A mountain tapir, puma and andean bear interacting model
R
occupancy
tapir
Author
Affiliation

Diego J. Lizcano

WildMon

Published

July 7, 2024

How can we model species interactions?:

  1. Direct observations of interactions (e.g. predation events)
  2. Indirect ways:

  • Over time. We can use 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.

  • Over space. Multispecies occupancy models. Don’t necessary test if species are active at the same time.

Three are several types of multispecies occupancy models:

  • Two or more 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 post.

Load packages

Code

# library(ggpmthemes)
library(glue) # Interpreted String Literals
library(patchwork) # The Composer of Plots
library(readxl) # Read Excel Files
library(sf) # Simple Features for R
library(mapview) # Interactive Viewing of Spatial Data in R
library(grateful) # Facilitate Citation of R Packages
library (terra)
library(unmarked)
library(stars)
library(elevatr)
library(ubms)
library(camtrapR)
library(knitr) # A General-Purpose Package for Dynamic Report Generation in R
# options(kableExtra.auto_format = FALSE)
library(kableExtra) # Construct Complex Table with 'kable' and Pipe Syntax
library(DT)
library(tidyverse) # Easily Install and Load the 'Tidyverse'
library(ggforce) # Accelerating 'ggplot2'

library(readr)



source("C:/CodigoR/CameraTrapCesar/R/organiza_datos.R")

Load data

Code

odp_ecu_col<- read_csv("C:/CodigoR/CameraTrapCesar/data/oso_danta_puma_ucu_pitalito_cocha2_ecuador_noNA.csv")

Convert to sf

Code

datos_distinct <- odp_ecu_col |> distinct(Longitude, Latitude, Deployment_id)

projlatlon <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

datos_sf <-  st_as_sf(x = datos_distinct,
                         coords = c("Longitude", 
                                    "Latitude"),
                         crs = projlatlon)

mapview(datos_sf, zcol="Deployment_id")

get rasters

Code

#load raster
per_tree_cov <- rast("F:/WCS-CameraTrap/raster/latlon/Veg_Cont_Fields_Yearly_250m_v61/Perc_TreeCov/MOD44B_Perc_TreeCov_2010_065.tif")
road_den <- rast("F:/WCS-CameraTrap/raster/latlon/RoadDensity/grip4_total_dens_m_km2.asc")
# elev <- rast("D:/CORREGIDAS/elevation_z7.tif")
landcov <- rast("F:/WCS-CameraTrap/raster/latlon/LandCover_Type_Yearly_500m_v61/LC1/MCD12Q1_LC1_2010_001.tif") 
cattle <- rast("F:/WCS-CameraTrap/raster/latlon/Global cattle distribution/5_Ct_2010_Da.tif")
#river <- st_read("F:/WCS-CameraTrap/shp/DensidadRios/MCD12Q1_LC1_2001_001_RECLASS_MASK_GRID_3600m_DensDrenSouthAmer.shp")

# get elevation map
#elevation_detailed <- rast(get_elev_raster(sites, z = 10, clip="bbox", neg_to_na=TRUE))
elevation_detailed <- get_elev_point (datos_sf, src="aws", overwrite=TRUE)


# extract covs using points and add to sites
# covs <- cbind(sites, terra::extract(SiteCovsRast, sites))
per_tre <- terra::extract(per_tree_cov, datos_sf)
roads <- terra::extract(road_den, datos_sf)
# eleva <- terra::extract(elevation_detailed, sites)
land_cov <- terra::extract(landcov, datos_sf)
cattle_den <-  terra::extract(cattle, datos_sf)

sites <- as.data.frame(datos_sf)

# remove decimals convert to factor
sites$land_cover <-  factor(land_cov$MCD12Q1_LC1_2010_001)
# sites$elevation <-  eleva$file3be898018c3
sites$per_tree_cov <- per_tre$MOD44B_Perc_TreeCov_2010_065 
#  fix 200 isue
ind <- which(sites$per_tree_cov== 200)
sites$per_tree_cov[ind] <- 0

sites$elevation <- elevation_detailed$elevation
sites$roads <- roads$grip4_total_dens_m_km2
sites$cattle <- cattle_den[,2]

# arrange detections observations

effort <- as.data.frame(odp_ecu_col[,19:23])

ylist <- list(oso  =as.matrix(odp_ecu_col[,4:8]),
              danta=as.matrix(odp_ecu_col[,9:13]),
              puma =as.matrix(odp_ecu_col[,14:18]))

lapply(ylist, head) # look at first few rows
#> $oso
#>      o1 o2 o3 o4 o5
#> [1,]  0  0  0  0  0
#> [2,]  0  0  0  0  0
#> [3,]  0  0  0  0  0
#> [4,]  0  0  0  0  0
#> [5,]  0  0  0  0  0
#> [6,]  0  0  0  0  0
#> 
#> $danta
#>      d1 d2 d3 d4 d5
#> [1,]  1  1  1  1  1
#> [2,]  1  1  1  1  1
#> [3,]  1  1  1  0  1
#> [4,]  1  1  1  1  1
#> [5,]  0  1  1  1  1
#> [6,]  1  1  1  1  1
#> 
#> $puma
#>      p1 p2 p3 p4 p5
#> [1,]  0  0  0  0  0
#> [2,]  0  0  0  0  0
#> [3,]  0  0  1  0  0
#> [4,]  0  0  0  0  0
#> [5,]  0  0  0  0  0
#> [6,]  0  1  0  0  0

site_covs <- as.data.frame(sites[,4:7])
# site_covs$roads <- as.numeric(site_covs$roads)
# site_covs <- scale(site_covs)
head(site_covs)
#>   per_tree_cov elevation roads   cattle
#> 1           19      2199   541 2554.303
#> 2           13      2729   541 2554.303
#> 3            9      2199   541 2554.303
#> 4            5      2199   541 2554.303
#> 5            8      2199   541 2554.303
#> 6           13      2370   541 2554.303


ObsCovs_list <- list(effort= odp_ecu_col[,19:23])


# Make UMF object
umf <- unmarkedFrameOccuMulti(y=ylist, 
                              siteCovs=site_covs,
                              obsCovs=ObsCovs_list
                              )
head(umf)
#> Data frame representation of unmarkedFrame object.
#> Only showing observation matrix for species 1.
#>    y.1 y.2 y.3 y.4 y.5 per_tree_cov elevation roads   cattle effort.1 effort.2
#> 1    0   0   0   0   0           19      2199   541 2554.303      9.5       10
#> 2    0   0   0   0   0           13      2729   541 2554.303      9.5       10
#> 3    0   0   0   0   0            9      2199   541 2554.303      9.5       10
#> 4    0   0   0   0   0            5      2199   541 2554.303      9.5       10
#> 5    0   0   0   0   0            8      2199   541 2554.303      9.5       10
#> 6    0   0   0   0   0           13      2370   541 2554.303      9.5       10
#> 7    0   0   0   0   0            4      2218    57 2545.696      9.5       10
#> 8    0   0   0   0   0            5      2310    57 2545.696      9.5       10
#> 9    0   0   0   0   0           10      2310    57 2545.696      9.5       10
#> 10   0   1   0   0   0            6      2023  1315 2963.785      9.5       10
#>    effort.3 effort.4 effort.5
#> 1        10       10      4.5
#> 2        10       10      2.5
#> 3        10       10      4.5
#> 4        10       10      4.5
#> 5        10       10     10.0
#> 6        10       10      3.5
#> 7        10       10     10.0
#> 8        10       10      2.5
#> 9        10       10     10.0
#> 10       10       10      2.5
summary(umf)
#> unmarkedFrame Object
#> 
#> 261 sites
#> 3 species: oso danta puma 
#> Maximum number of observations per site: 5 
#> Mean number of observations per site:
#> oso: 5  danta: 5  puma: 5  
#> Sites with at least one detection:
#> oso: 46  danta: 86  puma: 25  
#> Tabulation of y observations:
#> oso:
#>    0    1 
#> 1238   67 
#> danta:
#>    0    1 
#> 1111  194 
#> puma:
#>    0    1 
#> 1267   38 
#> 
#> Site-level covariates:
#>   per_tree_cov     elevation        roads            cattle      
#>  Min.   : 2.00   Min.   : 684   Min.   :   0.0   Min.   :   0.0  
#>  1st Qu.: 9.00   1st Qu.:2232   1st Qu.:   0.0   1st Qu.: 220.9  
#>  Median :19.00   Median :2590   Median :  59.0   Median :1106.0  
#>  Mean   :29.02   Mean   :2631   Mean   : 180.5   Mean   :1288.1  
#>  3rd Qu.:49.00   3rd Qu.:3077   3rd Qu.: 227.0   3rd Qu.:2545.7  
#>  Max.   :82.00   Max.   :4019   Max.   :1315.0   Max.   :5666.5  
#> 
#> Observation-level covariates:
#>      effort     
#>  Min.   : 0.00  
#>  1st Qu.: 9.50  
#>  Median :10.00  
#>  Mean   :10.58  
#>  3rd Qu.:10.50  
#>  Max.   :20.00

plot(umf)

Set up the formulas

Intercept-only model, assuming independence

For now, we assume independence among species. We do this by only allowing 1st order natural parameters (maxOrder = 1).

This is equivalent to fitting 3 single-species occupancy models.

Code

umf@fDesign
#>          f1[oso] f2[danta] f3[puma] f4[oso:danta] f5[oso:puma] f6[danta:puma]
#> psi[111]       1         1        1             1            1              1
#> psi[110]       1         1        0             1            0              0
#> psi[101]       1         0        1             0            1              0
#> psi[100]       1         0        0             0            0              0
#> psi[011]       0         1        1             0            0              1
#> psi[010]       0         1        0             0            0              0
#> psi[001]       0         0        1             0            0              0
#> psi[000]       0         0        0             0            0              0
#>          f7[oso:danta:puma]
#> psi[111]                  1
#> psi[110]                  0
#> psi[101]                  0
#> psi[100]                  0
#> psi[011]                  0
#> psi[010]                  0
#> psi[001]                  0
#> psi[000]                  0

fit_1 <- occuMulti(detformulas = c('~1', '~1', '~1'),
                   stateformulas = c('~1', '~1', '~1'),
                   maxOrder = 1,
                   data = umf)

# DT::datatable(round(summary(fit_1), 3))

summary(fit_1)
#> 
#> Call:
#> occuMulti(detformulas = c("~1", "~1", "~1"), stateformulas = c("~1", 
#>     "~1", "~1"), data = umf, maxOrder = 1)
#> 
#> Occupancy (logit-scale):
#>                     Estimate    SE     z  P(>|z|)
#> [oso] (Intercept)     -0.988 0.241 -4.10 4.08e-05
#> [danta] (Intercept)   -0.609 0.138 -4.40 1.08e-05
#> [puma] (Intercept)    -1.833 0.269 -6.82 9.09e-12
#> 
#> Detection (logit-scale):
#>                     Estimate    SE     z  P(>|z|)
#> [oso] (Intercept)     -1.455 0.223 -6.53 6.63e-11
#> [danta] (Intercept)   -0.314 0.109 -2.88 3.99e-03
#> [puma] (Intercept)    -1.318 0.285 -4.63 3.69e-06
#> 
#> AIC: 1725.777 
#> Number of sites: 261
#> optim convergence code: 0
#> optim iterations: 32 
#> Bootstrap iterations: 0

Intercept-only model, assuming dependence

  • Set maxOrder = 2 to estimate up to 2nd order natural parameters
  • Permits dependence between species
  • Fixes all natural parameters > maxOrder at 0
Code

fit_2 <- occuMulti(detformulas = c('~1', '~1', '~1'),
                   stateformulas = c('~1', '~1', '~1',
                                     '~1', '~1', '~1'),
                   maxOrder = 2,
                   data = umf)

summary(fit_2)
#> 
#> Call:
#> occuMulti(detformulas = c("~1", "~1", "~1"), stateformulas = c("~1", 
#>     "~1", "~1", "~1", "~1", "~1"), data = umf, maxOrder = 2)
#> 
#> Occupancy (logit-scale):
#>                          Estimate    SE       z  P(>|z|)
#> [oso] (Intercept)         -0.6163 0.288 -2.1435 3.21e-02
#> [danta] (Intercept)       -0.5882 0.202 -2.9101 3.61e-03
#> [puma] (Intercept)        -2.7369 0.524 -5.2255 1.74e-07
#> [oso:danta] (Intercept)   -1.2996 0.505 -2.5737 1.01e-02
#> [oso:puma] (Intercept)    -0.0574 0.792 -0.0724 9.42e-01
#> [danta:puma] (Intercept)   1.8097 0.574  3.1527 1.62e-03
#> 
#> Detection (logit-scale):
#>                     Estimate    SE     z  P(>|z|)
#> [oso] (Intercept)      -1.46 0.223 -6.52 6.92e-11
#> [danta] (Intercept)    -0.31 0.109 -2.86 4.27e-03
#> [puma] (Intercept)     -1.32 0.285 -4.63 3.71e-06
#> 
#> AIC: 1708.863 
#> Number of sites: 261
#> optim convergence code: 0
#> optim iterations: 38 
#> Bootstrap iterations: 0

oso y danta occur together more frequently than expected by chance (p<0.01) danta y puma occur together more frequently than expected by chance (p<0.01)

Incorporating covariates

Any parameter can be modeled as a function of covariates Covariate models for each parameter can be unique Names of detection covariates correspond to names provided in named list Names of occupancy covariates correspond to names in data.frame The model below is driven less by biology and more by an interest in demonstrating that each parameter can be modeled uniquely.

Code
fit_3 <- occuMulti(detformulas = c('~effort', '~effort', '~effort'),
                   stateformulas = c('~cattle', #oso
                                     '~per_tree_cov', #danta
                                     '~elevation', #puma
                                     '~per_tree_cov', #oso:danta
                                     '~1', #oso:puma
                                     '~per_tree_cov' #danta:puma
                                     ),
                   maxOrder = 2,
                   se=TRUE,
                   penalty=0.5,
                   data = umf)
#> Bootstraping covariance matrix

summary(fit_3)
#> 
#> Call:
#> occuMulti(detformulas = c("~effort", "~effort", "~effort"), stateformulas = c("~cattle", 
#>     "~per_tree_cov", "~elevation", "~per_tree_cov", "~1", "~per_tree_cov"), 
#>     data = umf, maxOrder = 2, penalty = 0.5, se = TRUE)
#> 
#> Occupancy (logit-scale):
#>                           Estimate       SE      z  P(>|z|)
#> [oso] (Intercept)          2.43291 1.391020  1.749 8.03e-02
#> [oso] cattle              -0.00199 0.000747 -2.658 7.85e-03
#> [danta] (Intercept)        0.48386 0.413108  1.171 2.41e-01
#> [danta] per_tree_cov      -0.02671 0.011359 -2.351 1.87e-02
#> [puma] (Intercept)         1.72385 0.649402  2.655 7.94e-03
#> [puma] elevation          -0.00174 0.000378 -4.590 4.42e-06
#> [oso:danta] (Intercept)   -0.11941 0.808593 -0.148 8.83e-01
#> [oso:danta] per_tree_cov  -0.07561 0.046987 -1.609 1.08e-01
#> [oso:puma] (Intercept)    -0.56580 0.636328 -0.889 3.74e-01
#> [danta:puma] (Intercept)   0.97265 0.698584  1.392 1.64e-01
#> [danta:puma] per_tree_cov  0.01900 0.014875  1.278 2.01e-01
#> 
#> Detection (logit-scale):
#>                     Estimate     SE     z  P(>|z|)
#> [oso] (Intercept)    -3.6925 1.0052 -3.67 0.000239
#> [oso] effort          0.1333 0.0567  2.35 0.018820
#> [danta] (Intercept)  -0.8891 0.3963 -2.24 0.024851
#> [danta] effort        0.0594 0.0379  1.57 0.117158
#> [puma] (Intercept)   -1.7578 0.6764 -2.60 0.009357
#> [puma] effort         0.0480 0.0464  1.03 0.301497
#> 
#> AIC: 1685.991 
#> Number of sites: 261
#> optim convergence code: 0
#> optim iterations: 128 
#> Bootstrap iterations: 30

Conditional occupancy probability

Calculation of conditional and marginal occupancy probabilities is done with the predict function.

Create a data.frame for predictions The procedure is equivalent to creating data frames for all other applications of predict Include complete range of observed cattle; hold all other variables at their mean.

Code
nd_cond <- data.frame(
  cattle = rep(mean(site_covs$cattle), 100),
  elevation = rep(mean(site_covs$elevation), 100),
  roads = rep(mean(site_covs$roads), 100),
  per_tree_cov = seq(min(site_covs$per_tree_cov), max(site_covs$per_tree_cov),
                 length.out = 100)
)

Predicting danta occurrence when oso are present

species indicates which species we assume when predicting occupancy cond indicates which species we are assuming is present or absent

Code
danta_oso_1 <- predict(fit_3, type = 'state', species = 'danta',
                     cond = 'oso', newdata = nd_cond)

Predicting danta occurrence when oso are absent

putting a - in front of oso tells predict you wish to assume oso are absent

Code

danta_oso_0 <- predict(fit_3, type = 'state', species = 'danta',
                     cond = '-oso', newdata = nd_cond)

danta oso marginal occupancy box plot

Code
################################## Marginal
danta_marginal <- predict(fit_3, type="state", species="danta")
head(danta_marginal)
#>   Predicted         SE     lower     upper
#> 1 0.5475922 0.06877945 0.4229338 0.6783526
#> 2 0.5490742 0.07095654 0.4157483 0.7011157
#> 3 0.6025791 0.07310459 0.4753362 0.7469048
#> 4 0.6253549 0.07867852 0.4951861 0.7790982
#> 5 0.6082291 0.07432900 0.4770212 0.7548960
#> 6 0.5679836 0.06953967 0.4420549 0.7134346

oso_marginal <- predict(fit_3, type='state', species="oso") # get coyote
marg_plot_dat <- rbind(danta_marginal[1,], oso_marginal[1,])
marg_plot_dat$Species <- c("Tapir", "Bear")
marg_plot_dat
#>    Predicted         SE       lower     upper Species
#> 1 0.54759219 0.06877945 0.422933786 0.6783526   Tapir
#> 2 0.03576262 0.03955719 0.008303715 0.1521894    Bear


plot(1:2, marg_plot_dat$Predicted, ylim=c(0,0.9), 
     xlim=c(0.5,2.5), pch=19, cex=1.5, xaxt='n', 
     xlab="", ylab="Marginal occupancy")
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)
}

danta puma marginal occupancy box plot

Code
################################## Marginal
danta_marginal <- predict(fit_3, type="state", species="danta")
head(danta_marginal)
#>   Predicted         SE     lower     upper
#> 1 0.5475922 0.07803240 0.3773334 0.6672646
#> 2 0.5490742 0.07715404 0.3829454 0.6721972
#> 3 0.6025791 0.08146052 0.4282122 0.7278124
#> 4 0.6253549 0.08546553 0.4366426 0.7528809
#> 5 0.6082291 0.08232832 0.4302772 0.7339994
#> 6 0.5679836 0.07793673 0.4021121 0.6874401

oso_marginal <- predict(fit_3, type='state', species="puma") # get coyote
marg_plot_dat <- rbind(danta_marginal[1,], oso_marginal[1,])
marg_plot_dat$Species <- c("Tapir", "Puma")
marg_plot_dat
#>   Predicted        SE     lower     upper Species
#> 1 0.5475922 0.0780324 0.3773334 0.6672646   Tapir
#> 2 0.2219948 0.0609948 0.1309248 0.3577091    Puma


plot(1:2, marg_plot_dat$Predicted, ylim=c(0,0.9), 
     xlim=c(0.5,2.5), pch=19, cex=1.5, xaxt='n', 
     xlab="", ylab="Marginal occupancy")
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)
}

danta oso conditional box plot

Code

######################### Conditional
danta_oso <- predict(fit_3, type="state", species="danta", cond="oso")
head(danta_oso)
#>   Predicted        SE      lower     upper
#> 1 0.1959403 0.1081286 0.07927054 0.4821200
#> 2 0.2884221 0.1246683 0.10577594 0.5619816
#> 3 0.3951732 0.1465305 0.15525014 0.6980794
#> 4 0.4926029 0.1685261 0.17261997 0.7768049
#> 5 0.4190464 0.1520591 0.15956233 0.7232008
#> 6 0.2985319 0.1257494 0.10944791 0.5712308

danta_No_oso <- predict(fit_3, type="state", species="danta", cond="-oso")
head(danta_oso)
#>   Predicted        SE      lower     upper
#> 1 0.1959403 0.1081286 0.07927054 0.4821200
#> 2 0.2884221 0.1246683 0.10577594 0.5619816
#> 3 0.3951732 0.1465305 0.15525014 0.6980794
#> 4 0.4926029 0.1685261 0.17261997 0.7768049
#> 5 0.4190464 0.1520591 0.15956233 0.7232008
#> 6 0.2985319 0.1257494 0.10944791 0.5712308


plot_data <- rbind(danta_oso[1,], danta_No_oso[1,])
plot_data$Oso_status <- c("Present","Absent")
head(plot_data)
#>   Predicted        SE      lower     upper Oso_status
#> 1 0.1959403 0.1081286 0.07927054 0.4821200    Present
#> 2 0.5606346 0.0833134 0.38075351 0.7185665     Absent

plot(1:2, plot_data$Predicted, ylim=c(0, 0.9), 
     xlim=c(0.5,2.5), pch=19, cex=1.5, xaxt='n', 
     xlab="Bear status", ylab="Tapir cond. occupancy")
axis(1, at=1:2, labels=plot_data$Oso_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)
}

danta puma conditional box plot

Code

######################### Conditional
danta_oso <- predict(fit_3, type="state", species="danta", cond="puma")
head(danta_oso)
#>   Predicted         SE     lower     upper
#> 1 0.7822699 0.09819067 0.5341260 0.9227781
#> 2 0.7908443 0.10068087 0.5368515 0.9277790
#> 3 0.7966546 0.10378875 0.5226327 0.9313509
#> 4 0.8025998 0.10792380 0.5278690 0.9359451
#> 5 0.7981259 0.10472926 0.5191097 0.9325056
#> 6 0.7908443 0.10068087 0.5368515 0.9277790

danta_No_oso <- predict(fit_3, type="state", species="danta", cond="-puma")
head(danta_oso)
#>   Predicted         SE     lower     upper
#> 1 0.7822699 0.09819067 0.5341260 0.9227781
#> 2 0.7908443 0.10068087 0.5368515 0.9277790
#> 3 0.7966546 0.10378875 0.5226327 0.9313509
#> 4 0.8025998 0.10792380 0.5278690 0.9359451
#> 5 0.7981259 0.10472926 0.5191097 0.9325056
#> 6 0.7908443 0.10068087 0.5368515 0.9277790


plot_data <- rbind(danta_oso[1,], danta_No_oso[1,])
plot_data$Oso_status <- c("Present","Absent")
head(plot_data)
#>   Predicted         SE     lower     upper Oso_status
#> 1 0.7822699 0.09819067 0.5341260 0.9227781    Present
#> 2 0.4806296 0.08118172 0.3123268 0.6274171     Absent

plot(1:2, plot_data$Predicted, ylim=c(0, 0.9), 
     xlim=c(0.5,2.5), pch=19, cex=1.5, xaxt='n', 
     xlab="Puma status", ylab="Tapir cond. occupancy")
axis(1, at=1:2, labels=plot_data$Oso_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)
}

predicting with covariates

Code

gg_df_cond1 <- data.frame(
  per_tree_cov = rep(nd_cond$per_tree_cov, 2),
  occupancy = c(danta_oso_1$Predicted,
                danta_oso_0$Predicted),
  low = c(danta_oso_1$lower,
          danta_oso_0$lower),
  high = c(danta_oso_1$upper,
           danta_oso_0$upper),
  conditional = rep(c('Bear present', 'Bear absent'),
                    each = 100)
)


cond_fig1 <- ggplot(gg_df_cond1, aes(x = per_tree_cov, y = occupancy,
                                   group = conditional)) +
  geom_ribbon(aes(ymin = low, ymax = high, fill = conditional),  alpha=0.5) +
  geom_line() +
  ylab('Conditional Tapir\noccupancy probability') +
  xlab('per_tree_cov') +
  labs(fill = 'Bear state') +
  theme(text = element_text(size = 15),
        legend.position = c(0.75, 0.85))

cond_fig1

danta puma

Code

danta_puma_1 <- predict(fit_3, type = 'state', species = 'danta',
                     cond = 'puma', newdata = nd_cond)

danta_puma_0 <- predict(fit_3, type = 'state', species = 'danta',
                     cond = '-puma', newdata = nd_cond)


gg_df_cond2 <- data.frame(
  per_tree_cov = rep(nd_cond$per_tree_cov, 2),
  occupancy = c(danta_puma_1$Predicted,
                danta_puma_0$Predicted),
  low = c(danta_puma_1$lower,
          danta_puma_0$lower),
  high = c(danta_puma_1$upper,
           danta_puma_0$upper),
  conditional = rep(c('Puma present', 'Puma absent'),
                    each = 100)
)



cond_fig2 <- ggplot(gg_df_cond2, aes(x = per_tree_cov, y = occupancy,
                                   group = conditional)) +
  geom_ribbon(aes(ymin = low, ymax = high, fill = conditional),  alpha=0.5) +
  geom_line() +
  ylab('Conditional Tapir\noccupancy probability') +
  xlab('per_tree_cov') +
  labs(fill = 'Puma state') +
  theme(text = element_text(size = 15),
        legend.position = c(0.75, 0.85))

cond_fig2

https://eesc.usgs.gov/MBR/workshops/ahm2023/04_Multispecies%20_occupancy/multispecies-occupancy.html

Package Citation

Code
pkgs <- cite_packages(output = "paragraph", out.dir = ".") #knitr::kable(pkgs)
pkgs

We used R version 4.3.2 (R Core Team 2023) and the following R packages: camtrapR v. 2.3.0 (Niedballa et al. 2016), devtools v. 2.4.5 (Wickham et al. 2022), DT v. 0.32 (Xie, Cheng, and Tan 2024), elevatr v. 0.99.0 (Hollister et al. 2023), ggforce v. 0.4.2 (Pedersen 2024a), glue v. 1.7.0 (Hester and Bryan 2024), kableExtra v. 1.4.0 (Zhu 2024), knitr v. 1.46 (Xie 2014, 2015, 2024), mapview v. 2.11.2 (Appelhans et al. 2023), patchwork v. 1.2.0 (Pedersen 2024b), quarto v. 1.4 (Allaire and Dervieux 2024), rmarkdown v. 2.27 (Xie, Allaire, and Grolemund 2018; Xie, Dervieux, and Riederer 2020; Allaire et al. 2024), sf v. 1.0.15 (Pebesma 2018; Pebesma and Bivand 2023a), stars v. 0.6.4 (Pebesma and Bivand 2023b), styler v. 1.10.3 (Müller and Walthert 2024), terra v. 1.7.71 (Hijmans 2024), tidyverse v. 2.0.0 (Wickham et al. 2019), ubms v. 1.2.6 (Kellner et al. 2021), unmarked v. 1.4.1 (Fiske and Chandler 2011; Kellner et al. 2023).

Sesion info

Session info 
#> ─ Session info ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
#>  setting  value
#>  version  R version 4.3.2 (2023-10-31 ucrt)
#>  os       Windows 10 x64 (build 19042)
#>  system   x86_64, mingw32
#>  ui       RTerm
#>  language (EN)
#>  collate  Spanish_Colombia.utf8
#>  ctype    Spanish_Colombia.utf8
#>  tz       America/Bogota
#>  date     2024-07-06
#>  pandoc   3.1.11 @ C:/Program Files/RStudio/resources/app/bin/quarto/bin/tools/ (via rmarkdown)
#> 
#> ─ Packages ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
#>  ! package           * version  date (UTC) lib source
#>    abind             * 1.4-5    2016-07-21 [1] CRAN (R 4.3.1)
#>    base64enc           0.1-3    2015-07-28 [1] CRAN (R 4.3.1)
#>    bit                 4.0.5    2022-11-15 [1] CRAN (R 4.3.2)
#>    bit64               4.0.5    2020-08-30 [1] CRAN (R 4.3.2)
#>    boot                1.3-28.1 2022-11-22 [2] CRAN (R 4.3.2)
#>    brew                1.0-10   2023-12-16 [1] CRAN (R 4.3.2)
#>    cachem              1.0.8    2023-05-01 [1] CRAN (R 4.3.2)
#>    camtrapR          * 2.3.0    2024-02-26 [1] CRAN (R 4.3.3)
#>    cellranger          1.1.0    2016-07-27 [1] CRAN (R 4.3.2)
#>    class               7.3-22   2023-05-03 [2] CRAN (R 4.3.2)
#>    classInt            0.4-10   2023-09-05 [1] CRAN (R 4.3.2)
#>    cli                 3.6.2    2023-12-11 [1] CRAN (R 4.3.2)
#>    codetools           0.2-19   2023-02-01 [2] CRAN (R 4.3.2)
#>    colorspace          2.1-0    2023-01-23 [1] CRAN (R 4.3.2)
#>    crayon              1.5.2    2022-09-29 [1] CRAN (R 4.3.2)
#>    crosstalk           1.2.1    2023-11-23 [1] CRAN (R 4.3.2)
#>    curl                5.2.0    2023-12-08 [1] CRAN (R 4.3.2)
#>    data.table          1.15.0   2024-01-30 [1] CRAN (R 4.3.2)
#>    DBI                 1.2.2    2024-02-16 [1] CRAN (R 4.3.2)
#>    devtools            2.4.5    2022-10-11 [1] CRAN (R 4.3.2)
#>    digest              0.6.34   2024-01-11 [1] CRAN (R 4.3.2)
#>    dplyr             * 1.1.4    2023-11-17 [1] CRAN (R 4.3.2)
#>    DT                * 0.32     2024-02-19 [1] CRAN (R 4.3.3)
#>    e1071               1.7-14   2023-12-06 [1] CRAN (R 4.3.2)
#>    elevatr           * 0.99.0   2023-09-12 [1] CRAN (R 4.3.2)
#>    ellipsis            0.3.2    2021-04-29 [1] CRAN (R 4.3.2)
#>    evaluate            0.23     2023-11-01 [1] CRAN (R 4.3.2)
#>    fansi               1.0.6    2023-12-08 [1] CRAN (R 4.3.2)
#>    farver              2.1.1    2022-07-06 [1] CRAN (R 4.3.2)
#>    fastmap             1.1.1    2023-02-24 [1] CRAN (R 4.3.2)
#>    forcats           * 1.0.0    2023-01-29 [1] CRAN (R 4.3.2)
#>    fs                  1.6.3    2023-07-20 [1] CRAN (R 4.3.2)
#>    generics            0.1.3    2022-07-05 [1] CRAN (R 4.3.2)
#>    ggforce           * 0.4.2    2024-02-19 [1] CRAN (R 4.3.3)
#>    ggplot2           * 3.5.1    2024-04-23 [1] CRAN (R 4.3.3)
#>    glue              * 1.7.0    2024-01-09 [1] CRAN (R 4.3.2)
#>    grateful          * 0.2.4    2023-10-22 [1] CRAN (R 4.3.3)
#>    gridExtra           2.3      2017-09-09 [1] CRAN (R 4.3.2)
#>    gtable              0.3.4    2023-08-21 [1] CRAN (R 4.3.2)
#>    hms                 1.1.3    2023-03-21 [1] CRAN (R 4.3.2)
#>    htmltools           0.5.7    2023-11-03 [1] CRAN (R 4.3.2)
#>    htmlwidgets         1.6.4    2023-12-06 [1] CRAN (R 4.3.2)
#>    httpuv              1.6.14   2024-01-26 [1] CRAN (R 4.3.2)
#>    httr                1.4.7    2023-08-15 [1] CRAN (R 4.3.2)
#>    inline              0.3.19   2021-05-31 [1] CRAN (R 4.3.2)
#>    jquerylib           0.1.4    2021-04-26 [1] CRAN (R 4.3.2)
#>    jsonlite            1.8.8    2023-12-04 [1] CRAN (R 4.3.2)
#>    kableExtra        * 1.4.0    2024-01-24 [1] CRAN (R 4.3.3)
#>    KernSmooth          2.23-22  2023-07-10 [2] CRAN (R 4.3.2)
#>    knitr             * 1.46     2024-04-06 [1] CRAN (R 4.3.3)
#>    labeling            0.4.3    2023-08-29 [1] CRAN (R 4.3.1)
#>    later               1.3.2    2023-12-06 [1] CRAN (R 4.3.2)
#>    lattice             0.22-5   2023-10-24 [1] CRAN (R 4.3.2)
#>    leafem              0.2.3    2023-09-17 [1] CRAN (R 4.3.2)
#>    leaflet             2.2.1    2023-11-13 [1] CRAN (R 4.3.2)
#>    leaflet.providers   2.0.0    2023-10-17 [1] CRAN (R 4.3.2)
#>    leafpop             0.1.0    2021-05-22 [1] CRAN (R 4.3.2)
#>    lifecycle           1.0.4    2023-11-07 [1] CRAN (R 4.3.2)
#>    lme4                1.1-35.3 2024-04-16 [1] CRAN (R 4.3.2)
#>    loo                 2.7.0    2024-02-24 [1] CRAN (R 4.3.2)
#>    lubridate         * 1.9.3    2023-09-27 [1] CRAN (R 4.3.2)
#>    magrittr            2.0.3    2022-03-30 [1] CRAN (R 4.3.2)
#>    mapview           * 2.11.2   2023-10-13 [1] CRAN (R 4.3.2)
#>    MASS                7.3-60   2023-05-04 [2] CRAN (R 4.3.2)
#>    Matrix              1.6-1.1  2023-09-18 [2] CRAN (R 4.3.2)
#>    matrixStats         1.2.0    2023-12-11 [1] CRAN (R 4.3.2)
#>    memoise             2.0.1    2021-11-26 [1] CRAN (R 4.3.2)
#>    mgcv                1.9-1    2023-12-21 [1] CRAN (R 4.3.3)
#>    mime                0.12     2021-09-28 [1] CRAN (R 4.3.1)
#>    miniUI              0.1.1.1  2018-05-18 [1] CRAN (R 4.3.2)
#>    minqa               1.2.6    2023-09-11 [1] CRAN (R 4.3.2)
#>    munsell             0.5.0    2018-06-12 [1] CRAN (R 4.3.2)
#>    nlme                3.1-163  2023-08-09 [2] CRAN (R 4.3.2)
#>    nloptr              2.0.3    2022-05-26 [1] CRAN (R 4.3.2)
#>    patchwork         * 1.2.0    2024-01-08 [1] CRAN (R 4.3.3)
#>    pbapply             1.7-2    2023-06-27 [1] CRAN (R 4.3.2)
#>    pillar              1.9.0    2023-03-22 [1] CRAN (R 4.3.2)
#>    pkgbuild            1.4.4    2024-03-17 [1] CRAN (R 4.3.3)
#>    pkgconfig           2.0.3    2019-09-22 [1] CRAN (R 4.3.2)
#>    pkgload             1.3.4    2024-01-16 [1] CRAN (R 4.3.2)
#>    png                 0.1-8    2022-11-29 [1] CRAN (R 4.3.1)
#>    polyclip            1.10-6   2023-09-27 [1] CRAN (R 4.3.1)
#>    prettyunits         1.2.0    2023-09-24 [1] CRAN (R 4.3.2)
#>    processx            3.8.3    2023-12-10 [1] CRAN (R 4.3.2)
#>    profvis             0.3.8    2023-05-02 [1] CRAN (R 4.3.2)
#>    progress            1.2.3    2023-12-06 [1] CRAN (R 4.3.3)
#>    progressr           0.14.0   2023-08-10 [1] CRAN (R 4.3.2)
#>    promises            1.2.1    2023-08-10 [1] CRAN (R 4.3.2)
#>    proxy               0.4-27   2022-06-09 [1] CRAN (R 4.3.2)
#>    ps                  1.7.6    2024-01-18 [1] CRAN (R 4.3.2)
#>    purrr             * 1.0.2    2023-08-10 [1] CRAN (R 4.3.2)
#>    quarto            * 1.4      2024-03-06 [1] CRAN (R 4.3.3)
#>    QuickJSR            1.1.3    2024-01-31 [1] CRAN (R 4.3.2)
#>    R.cache             0.16.0   2022-07-21 [1] CRAN (R 4.3.3)
#>    R.methodsS3         1.8.2    2022-06-13 [1] CRAN (R 4.3.3)
#>    R.oo                1.26.0   2024-01-24 [1] CRAN (R 4.3.3)
#>    R.utils             2.12.3   2023-11-18 [1] CRAN (R 4.3.3)
#>    R6                  2.5.1    2021-08-19 [1] CRAN (R 4.3.2)
#>    raster              3.6-26   2023-10-14 [1] CRAN (R 4.3.2)
#>    Rcpp                1.0.12   2024-01-09 [1] CRAN (R 4.3.2)
#>    RcppNumerical       0.6-0    2023-09-06 [1] CRAN (R 4.3.3)
#>  D RcppParallel        5.1.7    2023-02-27 [1] CRAN (R 4.3.2)
#>    readr             * 2.1.5    2024-01-10 [1] CRAN (R 4.3.2)
#>    readxl            * 1.4.3    2023-07-06 [1] CRAN (R 4.3.2)
#>    remotes             2.5.0    2024-03-17 [1] CRAN (R 4.3.3)
#>    renv                1.0.7    2024-04-11 [1] CRAN (R 4.3.3)
#>    rlang               1.1.3    2024-01-10 [1] CRAN (R 4.3.2)
#>    rmarkdown           2.27     2024-05-17 [1] CRAN (R 4.3.3)
#>    RSpectra            0.16-1   2022-04-24 [1] CRAN (R 4.3.2)
#>    rstan               2.32.6   2024-03-05 [1] CRAN (R 4.3.3)
#>    rstantools          2.4.0    2024-01-31 [1] CRAN (R 4.3.2)
#>    rstudioapi          0.16.0   2024-03-24 [1] CRAN (R 4.3.3)
#>    satellite           1.0.5    2024-02-10 [1] CRAN (R 4.3.2)
#>    scales              1.3.0    2023-11-28 [1] CRAN (R 4.3.3)
#>    secr                4.6.6    2024-02-29 [1] CRAN (R 4.3.3)
#>    sessioninfo         1.2.2    2021-12-06 [1] CRAN (R 4.3.2)
#>    sf                * 1.0-15   2023-12-18 [1] CRAN (R 4.3.2)
#>    shiny               1.8.0    2023-11-17 [1] CRAN (R 4.3.2)
#>    slippymath          0.3.1    2019-06-28 [1] CRAN (R 4.3.2)
#>    sp                  2.1-3    2024-01-30 [1] CRAN (R 4.3.2)
#>    StanHeaders         2.32.5   2024-01-10 [1] CRAN (R 4.3.2)
#>    stars             * 0.6-4    2023-09-11 [1] CRAN (R 4.3.2)
#>    stringi             1.8.3    2023-12-11 [1] CRAN (R 4.3.2)
#>    stringr           * 1.5.1    2023-11-14 [1] CRAN (R 4.3.2)
#>    styler            * 1.10.3   2024-04-07 [1] CRAN (R 4.3.3)
#>    svglite             2.1.3    2023-12-08 [1] CRAN (R 4.3.2)
#>    systemfonts         1.0.5    2023-10-09 [1] CRAN (R 4.3.2)
#>    terra             * 1.7-71   2024-01-31 [1] CRAN (R 4.3.2)
#>    tibble            * 3.2.1    2023-03-20 [1] CRAN (R 4.3.2)
#>    tidyr             * 1.3.1    2024-01-24 [1] CRAN (R 4.3.2)
#>    tidyselect          1.2.1    2024-03-11 [1] CRAN (R 4.3.3)
#>    tidyverse         * 2.0.0    2023-02-22 [1] CRAN (R 4.3.2)
#>    timechange          0.3.0    2024-01-18 [1] CRAN (R 4.3.2)
#>    tweenr              2.0.3    2024-02-26 [1] CRAN (R 4.3.3)
#>    tzdb                0.4.0    2023-05-12 [1] CRAN (R 4.3.2)
#>    ubms              * 1.2.6    2023-09-11 [1] CRAN (R 4.3.2)
#>    units               0.8-5    2023-11-28 [1] CRAN (R 4.3.2)
#>    unmarked          * 1.4.1    2024-01-09 [1] CRAN (R 4.3.2)
#>    urlchecker          1.0.1    2021-11-30 [1] CRAN (R 4.3.2)
#>    usethis             2.2.3    2024-02-19 [1] CRAN (R 4.3.2)
#>    utf8                1.2.4    2023-10-22 [1] CRAN (R 4.3.2)
#>    uuid                1.2-0    2024-01-14 [1] CRAN (R 4.3.2)
#>    V8                  4.4.2    2024-02-15 [1] CRAN (R 4.3.3)
#>    vctrs               0.6.5    2023-12-01 [1] CRAN (R 4.3.2)
#>    viridisLite         0.4.2    2023-05-02 [1] CRAN (R 4.3.2)
#>    vroom               1.6.5    2023-12-05 [1] CRAN (R 4.3.2)
#>    withr               3.0.0    2024-01-16 [1] CRAN (R 4.3.2)
#>    xfun                0.44     2024-05-15 [1] CRAN (R 4.3.3)
#>    xml2                1.3.6    2023-12-04 [1] CRAN (R 4.3.2)
#>    xtable              1.8-4    2019-04-21 [1] CRAN (R 4.3.2)
#>    yaml                2.3.8    2023-12-11 [1] CRAN (R 4.3.2)
#> 
#>  [1] C:/Users/usuario/AppData/Local/R/win-library/4.3
#>  [2] C:/Program Files/R/R-4.3.2/library
#> 
#>  D ── DLL MD5 mismatch, broken installation.
#> 
#> ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────

References

Allaire, JJ, and Christophe Dervieux. 2024. quarto: R Interface to “Quarto” Markdown Publishing System. https://CRAN.R-project.org/package=quarto.
Allaire, JJ, Yihui Xie, Christophe Dervieux, Jonathan McPherson, Javier Luraschi, Kevin Ushey, Aron Atkins, et al. 2024. rmarkdown: Dynamic Documents for r. https://github.com/rstudio/rmarkdown.
Appelhans, Tim, Florian Detsch, Christoph Reudenbach, and Stefan Woellauer. 2023. mapview: Interactive Viewing of Spatial Data in r. https://CRAN.R-project.org/package=mapview.
Fiske, Ian, and Richard Chandler. 2011. “unmarked: An R Package for Fitting Hierarchical Models of Wildlife Occurrence and Abundance.” Journal of Statistical Software 43 (10): 1–23. https://www.jstatsoft.org/v43/i10/.
Hester, Jim, and Jennifer Bryan. 2024. glue: Interpreted String Literals. https://CRAN.R-project.org/package=glue.
Hijmans, Robert J. 2024. terra: Spatial Data Analysis. https://CRAN.R-project.org/package=terra.
Hollister, Jeffrey, Tarak Shah, Jakub Nowosad, Alec L. Robitaille, Marcus W. Beck, and Mike Johnson. 2023. elevatr: Access Elevation Data from Various APIs. https://doi.org/10.5281/zenodo.8335450.
Kellner, Kenneth F., Nicholas L. Fowler, Tyler R. Petroelje, Todd M. Kautz, Dean E. Beyer, and Jerrold L. Belant. 2021. “ubms: An R Package for Fitting Hierarchical Occupancy and n-Mixture Abundance Models in a Bayesian Framework.” Methods in Ecology and Evolution 13: 577–84. https://doi.org/10.1111/2041-210X.13777.
Kellner, Kenneth F., Adam D. Smith, J. Andrew Royle, Marc Kery, Jerrold L. Belant, and Richard B. Chandler. 2023. “The unmarked R Package: Twelve Years of Advances in Occurrence and Abundance Modelling in Ecology.” Methods in Ecology and Evolution 14 (6): 1408–15. https://www.jstatsoft.org/v43/i10/.
Müller, Kirill, and Lorenz Walthert. 2024. styler: Non-Invasive Pretty Printing of r Code. https://CRAN.R-project.org/package=styler.
Niedballa, Jürgen, Rahel Sollmann, Alexandre Courtiol, and Andreas Wilting. 2016. “camtrapR: An r Package for Efficient Camera Trap Data Management.” Methods in Ecology and Evolution 7 (12): 1457–62. https://doi.org/10.1111/2041-210X.12600.
Pebesma, Edzer. 2018. “Simple Features for R: Standardized Support for Spatial Vector Data.” The R Journal 10 (1): 439–46. https://doi.org/10.32614/RJ-2018-009.
Pebesma, Edzer, and Roger Bivand. 2023a. Spatial Data Science: With applications in R. Chapman and Hall/CRC. https://doi.org/10.1201/9780429459016.
———. 2023b. Spatial Data Science: With applications in R. London: Chapman; Hall/CRC. https://doi.org/10.1201/9780429459016.
Pedersen, Thomas Lin. 2024a. ggforce: Accelerating “ggplot2”. https://CRAN.R-project.org/package=ggforce.
———. 2024b. patchwork: The Composer of Plots. https://CRAN.R-project.org/package=patchwork.
R Core Team. 2023. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.
Wickham, Hadley, Mara Averick, Jennifer Bryan, Winston Chang, Lucy D’Agostino McGowan, Romain François, Garrett Grolemund, et al. 2019. “Welcome to the tidyverse.” Journal of Open Source Software 4 (43): 1686. https://doi.org/10.21105/joss.01686.
Wickham, Hadley, Jim Hester, Winston Chang, and Jennifer Bryan. 2022. devtools: Tools to Make Developing r Packages Easier. https://CRAN.R-project.org/package=devtools.
Xie, Yihui. 2014. “knitr: A Comprehensive Tool for Reproducible Research in R.” In Implementing Reproducible Computational Research, edited by Victoria Stodden, Friedrich Leisch, and Roger D. Peng. Chapman; Hall/CRC.
———. 2015. Dynamic Documents with R and Knitr. 2nd ed. Boca Raton, Florida: Chapman; Hall/CRC. https://yihui.org/knitr/.
———. 2024. knitr: A General-Purpose Package for Dynamic Report Generation in r. https://yihui.org/knitr/.
Xie, Yihui, J. J. Allaire, and Garrett Grolemund. 2018. R Markdown: The Definitive Guide. Boca Raton, Florida: Chapman; Hall/CRC. https://bookdown.org/yihui/rmarkdown.
Xie, Yihui, Joe Cheng, and Xianying Tan. 2024. DT: A Wrapper of the JavaScript Library “DataTables”. https://CRAN.R-project.org/package=DT.
Xie, Yihui, Christophe Dervieux, and Emily Riederer. 2020. R Markdown Cookbook. Boca Raton, Florida: Chapman; Hall/CRC. https://bookdown.org/yihui/rmarkdown-cookbook.
Zhu, Hao. 2024. kableExtra: Construct Complex Table with “kable” and Pipe Syntax. https://CRAN.R-project.org/package=kableExtra.

Citation

BibTeX citation:
@online{j._lizcano2024,
  author = {J. Lizcano, Diego},
  title = {A Multi-Species (Species Interactions) Occupancy Model},
  date = {2024-07-07},
  url = {https://dlizcano.github.io/cametrapcesar/posts/2024-07-05-multi-species-interacting-occupancy/},
  langid = {en}
}
For attribution, please cite this work as:
J. Lizcano, Diego. 2024. “A Multi-Species (Species Interactions) Occupancy Model.” July 7, 2024. https://dlizcano.github.io/cametrapcesar/posts/2024-07-05-multi-species-interacting-occupancy/.