aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorKen Kellner <ken@kenkellner.com>2021-11-17 14:53:37 -0500
committerKen Kellner <ken@kenkellner.com>2021-11-17 14:53:37 -0500
commit3de220fb7d04b8825c46421c54ccc7117e03b7f8 (patch)
tree7578a22261c05344363e2b3f8da7968e14063167
parent319bb754a8ecdaf9470fbabf0902d6fc93af844b (diff)
Add occuMulti vignette
-rw-r--r--DESCRIPTION5
-rw-r--r--vignettes/occuMulti.Rnw518
-rw-r--r--vignettes/unmarked.bib19
3 files changed, 540 insertions, 2 deletions
diff --git a/DESCRIPTION b/DESCRIPTION
index 7bfb6b5..4ec8c96 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
Package: unmarked
-Version: 1.1.1.9013
-Date: 2021-11-08
+Version: 1.1.1.9014
+Date: 2021-11-17
Type: Package
Title: Models for Data from Unmarked Animals
Authors@R: c(
@@ -32,6 +32,7 @@ Imports:
stats,
TMB (>= 1.7.18),
utils
+Suggests: AHMbook
Description: Fits hierarchical models of animal abundance and occurrence to data collected using survey methods such as point counts, site occupancy sampling, distance sampling, removal sampling, and double observer sampling. Parameters governing the state and observation processes can be modeled as functions of covariates.
License: GPL (>=3)
LazyLoad: yes
diff --git a/vignettes/occuMulti.Rnw b/vignettes/occuMulti.Rnw
new file mode 100644
index 0000000..03ed0ad
--- /dev/null
+++ b/vignettes/occuMulti.Rnw
@@ -0,0 +1,518 @@
+<<echo=false>>=
+options(width=70)
+options(continue=" ")
+@
+
+\documentclass[a4paper]{article}
+\usepackage[OT1]{fontenc}
+\usepackage{Sweave}
+\usepackage{natbib}
+\usepackage{amsmath}
+%\usepackage{fullpage}
+\usepackage[vmargin=1in,hmargin=1in]{geometry}
+\bibliographystyle{ecology}
+
+\usepackage{hyperref}
+\hypersetup{
+ colorlinks=true,
+ linkcolor=blue,
+ urlcolor=cyan,
+ citecolor=black
+}
+
+\DefineVerbatimEnvironment{Sinput}{Verbatim} {xleftmargin=2em}
+\DefineVerbatimEnvironment{Soutput}{Verbatim}{xleftmargin=2em}
+\DefineVerbatimEnvironment{Scode}{Verbatim}{xleftmargin=2em}
+\fvset{listparameters={\setlength{\topsep}{0pt}}}
+\renewenvironment{Schunk}{\vspace{\topsep}}{\vspace{\topsep}}
+
+%%\VignetteIndexEntry{Multispecies occupancy models with occuMulti}
+
+\title{Multispecies occupancy models with occuMulti}
+\author{Ken Kellner}
+\date{November 15, 2021}
+
+
+\begin{document}
+
+\newcommand{\code}[1]{\texttt{\small{#1}}}
+\newcommand{\package}[1]{\textsf{\small{#1}}}
+
+\maketitle
+
+\section*{Outline}
+
+\begin{enumerate}
+ \item Introduction
+ \item Simple multispecies analysis
+ \item Analysis with covariates
+ \item Model selection
+ \item Model fitting challenges
+ \item Penalized likelihood
+\end{enumerate}
+
+
+\section{Introduction}
+
+The Rota et al. (\citeyear{Rota2016}) occupancy model is designed for presence/absence datasets with two or more (potentially) interacting species.
+The model allows for estimation of occupancy probabilities and the strength of interactions between species, as well as covariate effects on these parameters.
+The model generalizes the standard single-species occupancy model from MacKenzie et al. (\citeyear{mackenzie_estimating_2002}).
+
+The model assumes the latent occupancy state at site $i$ for a set of $s$ potentially interacting species is a vector $Z_i$ of length $s$ containing a sequence of the values 0 or 1.
+For example, when $s = 2$, the possible states are [11], [10], [01], or [00], corresponding to both species present, only species 1 or species 2 present, or both species absent, respectively.
+The latent state modeled as a multivariate Bernoulli random variable:
+\begin{equation}
+ Z_i \sim \mathrm{MVB}(\psi_i)
+\end{equation}
+where $\psi_i$ is a vector of length $2^s$ containing the probability of each possible combination of 0s and 1s, such that $\sum \psi_i = 1$.
+
+For $s = 2$, the corresponding natural parameters, $f$, are
+
+\begin{equation}
+\begin{split}
+ f_1 &= \mathrm{log}(\psi_{10}/\psi_{00}) \\
+ f_2 &= \mathrm{log}(\psi_{01}/\psi_{00}) \\
+ f_{12} &= \mathrm{log}((\psi_{11} * \psi_{00})/(\psi_{10} * \psi_{01}))
+\end{split}
+\end{equation}
+
+The natural parameters can then be modeled as linear functions of covariates.
+
+The observation process is similar to the standard single-species occupancy model, except that the observations $y_{ij}$ at site $i$ on occasion $j$ are vectors of length $s$ and there are independent values of detection probability $p$ for each species $s$:
+\begin{equation}
+ y_{ij} | Z_i \sim \mathrm{Bernoulli}(Z_i * p_{sij})
+\end{equation}
+See \citet{Rota2016} for more details on model structure.
+In \code{unmarked}, the model can be fit with the \code{occuMulti} function.
+
+\section{Simple multispecies analysis}
+
+We will start with a simple analysis using presence/absence data from camera traps for three species: coyote, red fox, and bobcat.
+The data are a simplified version of the data used in \citet{Rota2016}, with the data collapsed three into three-week long time periods.
+
+\subsection{Formatting the data}
+
+The dataset is called \code{MesoCarnivores} and is provided by the \code{AHMbook} package.
+First, we need to load in the dataset, which is a list with several components.
+
+<<>>=
+library(unmarked)
+library(AHMbook)
+data(MesoCarnivores)
+names(MesoCarnivores)
+@
+
+Presence/absence matrices for the three species are in list elements \code{bobcat}, \code{coyote}, and \code{redfox}, and \code{sitecovs} contains the site-level covariate data.
+
+Using this information, we will create an \code{unmarkedFrameOccuMulti} object.
+You can information by looking at the help file for \code{unmarkedFrameOccuMulti}:
+
+<<eval=FALSE>>=
+?unmarkedFrameOccuMulti
+@
+
+First we combine the detection data for the 3 species into one named list.
+These names will be used throughout the multispecies analysis, so pick useful ones.
+
+<<>>=
+ylist <- list(bobcat=MesoCarnivores$bobcat, coyote=MesoCarnivores$coyote,
+ redfox=MesoCarnivores$redfox)
+lapply(ylist, head)
+@
+
+The site covariates are contained in a data frame:
+
+<<>>=
+head(MesoCarnivores$sitecovs)
+@
+
+The site covariates include a metric of disturbance in a 5 km radius (\code{Dist\_5km}), housing density in a 5 km radius (\code{HDens\_5km}), latitude, longitude, a metric of how many people use the site (\code{People\_site}) and whether the camera site is on (or off) a trail.
+Using our \code{ylist}, site covariates, and observation covariates, we can construct an \code{unmarkedFrameOccuMulti} object.
+
+<<>>=
+umf <- unmarkedFrameOccuMulti(y=ylist, siteCovs=MesoCarnivores$sitecovs)
+@
+
+\subsection{Occupancy formulas}
+
+While most \code{unmarked} models have only one or two formulas (and a single formula for occupancy or abundance), \code{occuMulti} requires one formula per natural parameter $f$.
+Thus, there will be multiple formulas associated with occupancy.
+These formulas are organized into an ordered character vector.
+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 \code{unmarkedFrameOccuMulti}.
+
+<<>>=
+umf@fDesign
+@
+
+The number and order of the formulas in the vector should match the column names of this matrix.
+There are 7 columns in the matrix: thus, we'll need 7 formulas total, and they should be provided in the following order:
+
+<<>>=
+colnames(umf@fDesign)
+@
+
+For this model we'll set the 1st and 2nd-order $f$ parameters to be intercept-only, and fix the 3rd order parameter at 0.
+We will combine our formulas into a vector called \code{stateformulas}, like this:
+
+<<>>=
+stateformulas <- c("~1","~1","~1","~1","~1","~1","0")
+@
+
+Notice that the formulas are character strings (each wrapped in \code{""}). This is required.
+
+\subsection{Detection formulas}
+
+Each species has its own detection probability formula.
+Thus, there should be 3 total formulas combined in a \code{detformulas} vector.
+The order of the formulas should match the order of species in \code{ylist}.
+For this model, all three species will have intercept-only detection formulas.
+
+<<>>=
+detformulas <- c("~1","~1","~1")
+@
+
+\subsection{Fit the model}
+
+First, look at the help file for \code{occuMulti} to check what the required arguments are:
+
+<<eval=FALSE>>=
+?occuMulti
+@
+
+We now have all the pieces we need (\code{unmarkedFrameOccuMulti}, \code{stateformulas}, \code{detformulas}) needed to run a basic model which we will call \code{mod\_null}.
+
+<<>>=
+mod_null <- occuMulti(detformulas=detformulas, stateformulas=stateformulas, data=umf)
+summary(mod_null)
+@
+
+The regression parameters associated with each $f$ are identified by the species name (or combination of species names) in brackets.
+
+A few things to notice:
+
+\begin{itemize}
+ \item Coyote occupancy is the highest of the three species.
+ \item Negative relationship between bobcat and red fox
+ \item Positive relationships between coyote and the other two species
+ \item There is no three-species interaction term in the summary, because we fixed it at 0.
+\end{itemize}
+
+\subsection{Occupancy probabilities}
+
+To get the expected probability for each occupancy state ([11], [10] and so on) at each site, use the \code{predict} function.
+This gives you the probabilities along with standard errors and a 95% CI.
+
+<<>>=
+occ_prob <- predict(mod_null, type="state")
+head(occ_prob$Predicted)
+@
+
+The rows of this matrix should sum to 1.
+All rows are the same because estimated occupancies at all sites are the same - we didn't include any covariates.
+
+\subsection{Marginal occupancy}
+
+It's often more interesting to look at occupancy for species individually.
+For example, you might want to know the marginal occupancy of one species at each site (across all possible occupancy states).
+
+You can do this by specifying the \code{species} argument in \code{predict}, which will automatically sum up the appropriate occupancy states.
+
+<<>>=
+redfox_marginal <- predict(mod_null, type="state", species="redfox")
+head(redfox_marginal)
+@
+
+\subsection{Plotting marginal occupancy}
+
+Outputs from \code{predict} can be used to compare marginal occupancy across species with a plot.
+First, we'll need to get marginal occupancy for the other two species, and combine the three estimates into a single data frame.
+
+<<>>=
+coy_marginal <- predict(mod_null, type="state", species="coyote")
+bob_marginal <- predict(mod_null, type="state", species="bobcat")
+all_marginal <- rbind(redfox_marginal[1,], coy_marginal[1,], bob_marginal[1,])
+all_marginal$Species <- c("Red fox", "Coyote", "Bobcat")
+@
+
+Now we can plot the estimated marginal occupancy for each species, along with 95\% CIs.
+
+<<fig=TRUE>>=
+plot(1:3, all_marginal$Predicted, ylim=c(0.1,0.4),
+ xlim=c(0.5,3.5), pch=19, cex=1.5, xaxt='n',
+ xlab="", ylab="Marginal occupancy and 95% CI")
+axis(1, at=1:3, labels=all_marginal$Species)
+
+# CIs
+top <- 0.1
+for (i in 1:3){
+ segments(i, all_marginal$lower[i], i, all_marginal$upper[i])
+ segments(i-top, all_marginal$lower[i], i+top)
+ segments(i-top, all_marginal$upper[i], i+top)
+}
+@
+
+\subsection{Conditional occupancy}
+
+Alternatively, you might want to know the probability of occupancy of one species, conditional on the presence of another. Use the \code{species} and \code{cond} arguments together for this.
+For example, the probability of red fox occupancy, conditional on coyote presence:
+
+<<>>=
+redfox_coy <- predict(mod_null, type="state", species="redfox", cond="coyote")
+head(redfox_coy)
+@
+
+What about conditional on coyote *absence*?
+Simply add a \code{"-"} before the species name.
+
+<<>>=
+redfox_nocoy <- predict(mod_null, type="state", species="redfox", cond="-coyote")
+head(redfox_nocoy)
+@
+
+\subsection{Plotting conditional occupancy}
+
+You can use this output from \code{predict} to generate comparison plots.
+
+<<fig=TRUE>>=
+cond_data <- rbind(redfox_coy[1,], redfox_nocoy[1,])
+cond_data$Coyote_status <- c("Present","Absent")
+
+plot(1:2, cond_data$Predicted, ylim=c(0,0.3),
+ xlim=c(0.5,2.5), pch=19, cex=1.5, xaxt='n',
+ xlab="Coyote status", ylab="Red fox occupancy and 95% CI")
+axis(1, at=1:2, labels=cond_data$Coyote_status)
+
+# CIs
+top <- 0.1
+for (i in 1:2){
+ segments(i, cond_data$lower[i], i, cond_data$upper[i])
+ segments(i-top, cond_data$lower[i], i+top)
+ segments(i-top, cond_data$upper[i], i+top)
+}
+@
+
+Note that red fox occupancy is higher at sites where coyotes were present, which corresponds with the positive interaction term between the two species we saw in the model output summary.
+
+\section{Multispecies model with covariates}
+
+Now we'll fit a model with covariates on some natural parameters.
+It might be helpful to look at the order for our $f$ parameters again:
+
+<<>>=
+colnames(umf@fDesign)
+@
+
+And our available site covariates:
+
+<<>>=
+head(siteCovs(umf))
+@
+
+\subsection{Add housing density as a covariate}
+
+We'll fit a model with an effect of housing density (\code{HDens\_5km}) on the first-order parameters for all three species.
+The two-way interactions will remain intercept-only, and the three-way interaction will remain fixed at 0.
+Here's the vector of $f$ formulas:
+
+<<>>=
+sf <- c("~HDens_5km","~HDens_5km","~HDens_5km","~1","~1","~1","0")
+@
+
+Inside your formula, you can wrap the variable name in \code{scale()} to standardize it, which we should do, because the housing density variable is not scaled.
+
+<<>>=
+sf <- c("~scale(HDens_5km)","~scale(HDens_5km)","~scale(HDens_5km)","~1","~1","~1","0")
+@
+
+Detection formulas will remain the same, so we're now ready to fit a new model, \code{mod\_hdens}.
+
+<<>>=
+mod_hdens <- occuMulti(stateformulas=sf, detformulas=detformulas, umf)
+summary(mod_hdens)
+@
+
+A few things to note from the results:
+
+\begin{itemize}
+ \item Housing density has a significant negative effect on occupancy of bobcat
+ \item Housing density has a significant positive effect on red fox
+ \item No effect of housing density on coyote.
+\end{itemize}
+
+\subsection{Plotting covariate effects}
+
+To plot the effect of housing density on marginal occupancy, we again use \code{predict}.
+First, we need to generate sequence of possible \code{Hdens\_5km} values for the X-axis of our plot.
+We'll generate a sequence of 100 values, starting at the minimum observed housing density and ending at the maximum observed value.
+
+<<>>=
+hdens_range <- range(siteCovs(umf)$HDens_5km)
+hdens_seq <- seq(hdens_range[1], hdens_range[2], length.out=100)
+@
+
+Next, we'll \code{predict} marginal coyote occupancy at each value of \code{Hdens\_5km} along our sequence.
+Our custom housing density values should be in a data frame and provided to the \code{newdata} argument.
+Because we used \code{scale()} in our formulas above, there is no need to manually scale these new housing density values - \code{unmarked} will do it for us.
+
+<<>>=
+nd <- data.frame(HDens_5km = hdens_seq)
+occ_hdens_coy <- predict(mod_hdens, type="state", species="coyote", newdata=nd)
+occ_hdens_coy$Species <- "Coyote"
+occ_hdens_coy$Hdens <- hdens_seq
+head(occ_hdens_coy)
+@
+
+We'll do the same thing for the other two species.
+
+<<>>=
+occ_hdens_bob <- predict(mod_hdens, type="state", species="bobcat", newdata=nd)
+occ_hdens_fox <- predict(mod_hdens, type="state", species="redfox", newdata=nd)
+occ_hdens_bob$Species <- "Bobcat"
+occ_hdens_fox$Species <- "Red fox"
+occ_hdens_bob$Hdens <- hdens_seq
+occ_hdens_fox$Hdens <- hdens_seq
+@
+
+Finally, we'll build our plot.
+Housing density will be on the x-axis, marginal occupancy on the y-axis, and species will be identified by colors.
+
+<<fig=TRUE>>=
+plot(occ_hdens_coy$Hdens, occ_hdens_coy$Predicted, type='l', ylim=c(0,0.6),
+ col='red', lwd=2, xlab="Housing density", ylab="Marginal occupancy")
+lines(occ_hdens_bob$Hdens, occ_hdens_bob$Predicted, col='blue', lwd=2)
+lines(occ_hdens_fox$Hdens, occ_hdens_fox$Predicted, col='green', lwd=2)
+legend('topleft', col=c('red', 'blue', 'green'), lty=1,
+ legend=c("Coyote", "Bobcat", "Red fox"))
+@
+
+\section{Model selection}
+
+\code{unmarked} can calculate AIC, $\Delta$AIC, and weights for a list of models automatically.
+Start by creating a \code{fitList} object containing our two models:
+
+<<>>=
+mods <- fitList(mod_null, mod_hdens)
+@
+
+Then call the function \code{modSel} on our list of models to generate a model selection table:
+
+<<>>=
+modSel(mods)
+@
+
+Looks like the model with housing density is better.
+
+\section{Model fitting challenges}
+
+Multispecies occupancy models often have many parameters, and can be difficult to fit in some situations.
+You might get poor estimates (i.e., very large absolute values and/or large SEs) under certain conditions:
+
+\begin{itemize}
+ \item Sparse data (many 0s)
+ \item Boundary estimates (occupancy close to 0 or 1)
+ \item Few observations where multiple species are detected
+ \item Separation (perfect correlation with covariate)
+\end{itemize}
+
+Here's an example of a complex model with many covariates that results in poor estimates.
+
+<<>>=
+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)
+summary(mod_complex)
+@
+
+Note that several estimates are very large (>10) and also have large SEs.
+You should be very skeptical about using a model with poor estimates, like this one, for inference.
+Potential solutions when you get poor estimates include:
+
+\begin{itemize}
+ \item Fit simpler models with fewer covariates
+ \item If possible, fit the model with fewer species
+ \item Adjust observation period length if possible (e.g. collapse from one-week periods to three-week periods)
+ \item Use penalized likelihood to fit the model
+\end{itemize}
+
+\section{Penalized likelihood}
+
+\code{occuMulti} uses maximum likelihood to estimate parameters.
+We can add a "penalty" to the calculated likelihood to keep parameter estimates from getting stuck at huge values
+Use of penalized likelihiood has been shown to help with separation/boundary issues, eliminate unreasonably large estimates, and reduce error,
+However, note that the penalty term introduces a small amount of bias in the parameter estimates: thus, we are making a tradeoff between bias and variance
+With huge SEs as with the previous model, this may be a good tradeoff to make.
+
+One type of penalty is the "Bayes" penalty:
+\begin{equation}
+-\lambda\frac{1}{2}\sum_i{}\theta_i^2
+\end{equation}
+In the formula above, $\lambda$ is the penalty value, and $\theta$ is the vector of estimated parameters.
+As the parameter values get bigger, the total penalty increases.
+
+\subsection{Penalized likelihood with occuMulti}
+
+\code{occuMulti} can use penalized likelihood to fit models.
+You can provide a value to the \code{penalty} argument directly, or use the \code{optimizePenalty} function on a fitted model to choose the best value of $\lambda$ using K-fold cross-validation, and re-fit the model using the optimal penalty term value.
+
+<<eval=FALSE>>=
+set.seed(123)
+mod_penalty <- optimizePenalty(mod_complex, penalties=c(0.5,1))
+summary(mod_penalty)
+@
+
+<<eval=FALSE>>=
+## Optimal penalty is 1
+## Bootstraping covariance matrix
+
+## Call:
+## occuMulti(detformulas = c("~Trail", "~Trail", "~Trail"), stateformulas = c("~scale(Dist_5km)+scale(HDens_5km)",
+## "~scale(Dist_5km)+scale(HDens_5km)", "~scale(Dist_5km)+scale(HDens_5km)",
+## "~scale(Dist_5km)+scale(HDens_5km)", "~scale(Dist_5km)+scale(HDens_5km)",
+## "~scale(Dist_5km)+scale(HDens_5km)", "0"), data = object@data,
+## maxOrder = 3L, penalty = 1, boot = boot)
+
+## Occupancy (logit-scale):
+## Estimate SE z P(>|z|)
+## [bobcat] (Intercept) -1.7810 0.221 -8.054 8.00e-16
+## [bobcat] scale(Dist_5km) -1.3143 0.337 -3.903 9.48e-05
+## [bobcat] scale(HDens_5km) -2.8200 0.529 -5.334 9.61e-08
+## [coyote] (Intercept) -0.6049 0.178 -3.407 6.56e-04
+## [coyote] scale(Dist_5km) 0.0285 0.150 0.190 8.49e-01
+## [coyote] scale(HDens_5km) -1.0908 0.397 -2.748 5.99e-03
+## [redfox] (Intercept) -1.5659 0.310 -5.059 4.22e-07
+## [redfox] scale(Dist_5km) -0.3068 0.138 -2.226 2.60e-02
+## [redfox] scale(HDens_5km) 0.4730 0.797 0.593 5.53e-01
+## [bobcat:coyote] (Intercept) 1.1871 0.372 3.195 1.40e-03
+## [bobcat:coyote] scale(Dist_5km) 0.9347 0.368 2.537 1.12e-02
+## [bobcat:coyote] scale(HDens_5km) -0.3218 1.043 -0.309 7.58e-01
+## [bobcat:redfox] (Intercept) -0.8831 0.346 -2.553 1.07e-02
+## [bobcat:redfox] scale(Dist_5km) 0.0364 0.233 0.156 8.76e-01
+## [bobcat:redfox] scale(HDens_5km) 2.5609 1.074 2.384 1.71e-02
+## [coyote:redfox] (Intercept) 1.0001 0.249 4.009 6.09e-05
+## [coyote:redfox] scale(Dist_5km) 0.0236 0.229 0.103 9.18e-01
+## [coyote:redfox] scale(HDens_5km) 1.3920 0.424 3.281 1.03e-03
+
+## Detection (logit-scale):
+## Estimate SE z P(>|z|)
+## [bobcat] (Intercept) -2.44 0.150 -16.22 3.72e-59
+## [bobcat] Trail 1.74 0.164 10.61 2.59e-26
+## [coyote] (Intercept) -1.89 0.117 -16.24 2.72e-59
+## [coyote] Trail 2.10 0.150 14.00 1.52e-44
+## [redfox] (Intercept) -1.49 0.206 -7.21 5.66e-13
+## [redfox] Trail 1.72 0.254 6.79 1.14e-11
+
+## AIC: 6135.555
+## Number of sites: 1437
+## optim convergence code: 0
+## optim iterations: 100
+## Bootstrap iterations: 30
+@
+
+Notice that parameter estimates and SEs for the model using penalized likelihood are no longer gigantic.
+For more info on the use of penalized likelihood with multispecies occupancy models, see \citet{Clipp_2021}.
+
+\bibliography{unmarked}
+
+\end{document}
diff --git a/vignettes/unmarked.bib b/vignettes/unmarked.bib
index 672a401..b32d258 100644
--- a/vignettes/unmarked.bib
+++ b/vignettes/unmarked.bib
@@ -347,3 +347,22 @@ year = {2012}
}
+@article{Rota2016,
+ author = {Rota, Christopher T. and Ferreira, Marco A. R. and Kays, Roland and Forrester, Tavis D. and Kalies, Elizabeth L. and McShea, William J. and Parsons, Arielle W. and Millspaugh, Joshua J.},
+ title = {A multi-species occupancy model for two or more interacting species},
+ journal = {Methods in Ecology and Evolution},
+ year = {2016},
+ url = {http://dx.doi.org/10.1111/2041-210X.12587},
+ doi = {10.1111/2041-210X.12587},
+ volume = {7},
+ pages = {1164-1173},
+}
+
+@article{Clipp_2021,
+ author = {Clipp, Hannah L. and Evans, Amber L. and Kessinger, Brin E. and Kellner, Kenneth and Rota, Christopher T.},
+ title = {A penalized likelihood for multispecies occupancy models improves predictions of species interactions},
+ journal = {Ecology},
+ year = {2021},
+ pages = {e03520},
+ addendum = {\href{https://doi.org/10.1002/ecy.3520}{[view]}}
+}