# Assessing species interactions using integrated predator-prey models

Matthieu Paquet, Frederic Barraquand

► To cite this version:

Matthieu Paquet, Frederic Barraquand. Assessing species interactions using integrated predator-prey models. Peer Community Journal, 2023, 3, pp.e105. 10.24072/pcjournal.337. hal-04293759

HAL Id: hal-04293759

https://hal.science/hal-04293759v1

Submitted on 5 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
RESEARCH ARTICLE

Published

2023-11-08

Cite as

Matthieu Paquet and Frédéric Barraquand (2023) Assessing species interactions using integrated predator-prey models, Peer Community Journal, 3: e105.

Correspondence

matthieu.paquet@sete.cnrs.fr
freric.barraquand@u-bordeaux.fr

Peer-review

Peer reviewed and
recommended by
PCI Ecology,

https://doi.org/10.24072/pci.ecology.100522

(cc)

BY

This article is licensed
under the Creative Commons
Attribution 4.0 License.

Assessing species interactions
using integrated predator-prey
models

Matthieu Paquet<sup>1,0</sup> and Frédéric Barraquand<sup>1,0</sup>

Volume 3 (2023), article e105

https://doi.org/10.24072/pcjournal.337

Abstract

Inferring the strength of species interactions from demographic data is a challenging task. The Integrated Population Modelling (IPM) approach, bringing together population counts, capture-recapture, and individual-level fecundity data into a unified model framework, has been extended from single species to the community level. This allows to specify IPMs for multiple species with interactions specified as links between vital rates and stage-specific densities. However, there is no evaluation of such models when interactions are actually absent—while any interaction inference method runs the risk of producing false positives. We investigate here whether multispecies IPMs could output interactions where there are in fact none, building on an existing predator-prey IPM. We show that interspecific density-dependence estimates are centered on zero when simulated to be zero, and therefore their estimation is unbiased. Their coverage probability, quantifying how many times credible intervals include zero, is also satisfactory. We further confirm that adding random temporal variation to multispecies density-dependent link functions does not alter these results. This study therefore reaffirms the potential of multispecies IPMs to infer correctly how biotic interactions influence demography, although future studies should investigate model misspecifications.

<figure><img src="image_2.png" /></figure>

<footnote n="1">¹Institute of Mathematics of Bordeaux, University of Bordeaux, CNRS, Bordeaux INP, Talence, France</footnote>
Contents

<table><tr><td>1</td><td>Introduction</td><td>2</td></tr><tr><td>2</td><td>Methods</td><td>3</td></tr><tr><td>2.1</td><td>General description of the multispecies IPM</td><td>3</td></tr><tr><td>2.2</td><td>Alternative scenarios and parameter values</td><td>6</td></tr><tr><td>2.3</td><td>Prior specification and model fitting</td><td>8</td></tr><tr><td>3</td><td>Results</td><td>8</td></tr><tr><td>4</td><td>Discussion</td><td>13</td></tr><tr><td></td><td>Acknowledgements</td><td>15</td></tr><tr><td></td><td>Data and code</td><td>15</td></tr><tr><td></td><td>Funding</td><td>15</td></tr><tr><td></td><td>Conflict of interest disclosure</td><td>15</td></tr><tr><td></td><td>References</td><td>15</td></tr><tr><td></td><td>Supplementary Information</td><td>18</td></tr><tr><td>A</td><td>Results for the scenarios with species interactions</td><td>18</td></tr><tr><td>B</td><td>Results for the scenarios with 100 juveniles of each species marked each year for 10 years, and 20 juveniles of each species marked for 30 years, without centering abundances in the link functions</td><td>22</td></tr><tr><td>C</td><td>Sensitivity of parameter estimation to the choice of initial values</td><td>27</td></tr></table>

1. Introduction

Estimating ecological interactions between and within species through models of their joint population dynamics is a task which requires large amounts of data. Indeed, with potentially as many as $q^2$ interaction parameters for $q$ model compartments (combination of species and age classes), the number of parameters to estimate can climb very rapidly. Therefore, ecological statistics searches for improved ways to infer such population-level interaction strengths. A recently developed technique consists in combining data sources in multispecies Integrated Population Models (IPMs) including interspecific interactions (Barraquand and Gimenez, 2019; Péron and Koons, 2012; Quéroué et al., 2021). Because Integrated Population Models (IPMs, Besbeas et al., 2002) combine data on demographic rates (e.g., capture-recapture, breeding data) with data on population size (typically from counts), they allow: (a) estimating both demographic rates and population size (and hence their inter-dependencies) in a joint analysis, (b) an improved precision of parameter estimates, compared to separate analyses of component datasets, since the information contained in several datasets combine into estimated parameters (e.g., count data and capture recapture data both contain information on survival rates), and in some cases (c) to estimate parameters for which there is no dedicated data stream, that can only be estimated through inverse estimation of a demographic model (Abadi et al., 2010; Kéry and Schaub, 2011). This last property is particularly useful to estimate population-level species interaction strengths, since population-level interactions are always indirectly inferred. Although inverse estimation can in theory be performed using a single data source such as population counts, such inverse estimation is a difficult task fraught with identifiability issues. Asking whether multispecies IPMs performed better than classical inverse estimation from count data alone, Barraquand and Gimenez (2019) have shown that better estimates of interaction parameters could be obtained by combining data sources. Additionally, an empirical study in a bird predator-prey system (Quéroué et al., 2021) was able to detect the expected bottom-up demographic linkages from prey to predator but not the expected top-down relationships, suggesting that those may be too weak to be detected.

In these multispecies IPM studies estimating interspecific interactions, between-species linkages have always been considered to be present in the simulations or in the underlying reality (based on background knowledge). Other choices are possible: some multispecies IPMs do not
assume interspecific interactions to be present a priori (Lahoz-Monfort et al., 2017), but they do not estimate them either and focus instead on environmental effects. However, multispecies IPMs with interspecific interactions could also be used in situations where it is not clear whether population-level interactions between species are possible. This is all the more true that interactions are specified as links between vital rates and stage-specific densities, and while some of these relationships may be known a priori, others may not. The issue was raised but not tackled by Barraquand and Gimenez (2019): a natural follow-up is therefore to ask what happens whenever we try to estimate interactions that are actually absent, to make sure that multispecies IPMs do not yield false positives.

Let us note that when estimating or predicting interspecific interactions in general—not just with multispecies IPMs—whether methods could output false positives is a key concern (e.g., with multivariate autoregressive models, Barraquand et al. 2021; Mutshinda et al. 2009; dynamic bayesian networks, Sander et al. 2017; or other machine learning tools, Strydom et al. 2021). The fact that all interaction inference methods run the risk of creating false positives of interspecific interactions at exaggerated rates only reinforces the need to evaluate it in multispecies IPMs.

An additional concern is temporal stochasticity in the functions linking vital rates of a given stage of species *i* to the densities of a given stage of species *j*. In the simulation-based study of Barraquand and Gimenez (2019), it was assumed that such stochasticity was absent, while empirical studies (Péron and Koons, 2012; Quéroué et al., 2021) assumed its presence in order to partition variation in vital rates due to species densities vs other factors changing over time. We therefore still need to understand whether theoretical performances hold in this more empirically realistic context, where environmental factors can perturb demographic rates, and those are not solely deterministic functions of species densities.

To sum up, we follow-up here on the multispecies IPM study of Barraquand and Gimenez (2019) by asking whether (1) inter-species interactions are truly estimated to be zero when species have in fact independent dynamics and (2) how species interaction strengths estimates can be affected by the absence and presence of environmental stochasticity (random year effects on demographic rates).

## 2. Methods

### 2.1. General description of the multispecies IPM

The deterministic skeleton can be described as a density-dependent matrix population model (1)

$$ \mathbf{n}_{t+1} = \mathbf{A}(\mathbf{n}_t)\mathbf{n}_t. $$

Eq. 1 describes in discrete-time the dynamics of abundances of two species and two stages per species, with projection matrix

$$ \mathbf{A}(\mathbf{n}_t) = \begin{pmatrix} 0 & \frac{1}{2} f_{V,t} (n_{V,t}^A) \phi_{V,t}^J (n_{P,t}^A) & 0 & 0 \\ \phi_{V,t}^A & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{1}{2} f_{P,t} (n_{V,t}^J) \phi_{P,t}^J (n_{P,t}^A) \\ 0 & 0 & \phi_{P,t}^A & 0 \end{pmatrix} $$

and abundance vector

$$ \mathbf{n}_t = \begin{pmatrix} n_{V,t}^J \\ n_{V,t}^A \\ n_{P,t}^J \\ n_{P,t}^A \end{pmatrix} $$

where $n_{V,t}^J$, $n_{V,t}^A$, $n_{P,t}^J$ and $n_{P,t}^A$ are respectively the abundances of juvenile prey (denoted $V$ as ‘victim’), adult prey, juvenile predators and adult predators, at time $t$. The fecundities $f_{V,t}$, $f_{P,t}$ are the expected number of juvenile prey and predator produced by an adult female prey and predator, respectively. Survival probabilities between $t$ and $t+1$ are denoted with $\phi$, so that $\phi_{V,t}^J$, $\phi_{V,t}^A$, $\phi_{P,t}^J$ and $\phi_{P,t}^A$ are the survival probabilities of the juvenile prey, adult prey, juvenile predator and adult predator.
2.1.1. Count data. To simulate and account for demographic stochasticity, we modelled yearly (st)age specific abundances $n_t$ using Binomial and Poisson distributions as in Barraquand and Gimenez (2019) eqs. (2)–(5).

Regarding the observation process for count data, the 2019 model assumed a negligible observation error ($\sigma^2 = 10^{-5}$). The reason was that in absence of replicated counts at each time unit, observation error variance is notoriously difficult to disentangle from process error variance (Auger-Méthé et al., 2016; Knape, 2008). While in some cases it could be possible to remove observation error altogether, because total population sizes of each species (summed numbers of juveniles and adults) are the observed count variables (as in most IPMs), they need to appear in the model as drawn from some probability distribution—they need to be a stochastic node in the MCMC representation. It was therefore decided to keep the formulation of the model in its state-space version, but forcing it to observe true population size almost with certainty (negligible process error variance). However, we uncovered in the present work that stage-specific abundances could not be estimated properly. Because correctly reproducing stage-specific abundances when fitting a stage-structured model is desirable, and that there is in most wildlife surveys some measure of observation error on counts, we assumed in the present article a non-negligible, positive observation error variance. As we do not have replicated counts at any given time, we do not attempt to estimate observation error variance, and assume that it is known and classically set on the logarithmic scale (i.e., the coefficient of variation of observed population size is constant). For predator counts (denoted $P$) we have

$$ (2) \qquad y_{P,t} | \mathbf{n}_t \sim \mathcal{LN} \left( \log(n_{P,t}^J + n_{P,t}^A), \sigma_{obs}^2 \right) $$

and similarly for prey counts

$$ (3) \qquad y_{V,t} | \mathbf{n}_t \sim \mathcal{LN} \left( \log(n_{V,t}^J + n_{V,t}^A), \sigma_{obs}^2 \right), $$

with $\mathcal{LN}$ the log-Normal distribution and its associated standard deviation on the log-scale $\sigma_{obs} = 0.1$. Because $CV_{obs} = \sqrt{e^{\sigma_{obs}^2} - 1}$, this corresponds to $CV_{obs} \approx 10\%$. Other choices of observation model are possible but this one is standard for abundance values that are not too small (Besbeas et al., 2002; Dennis et al., 2006).

2.1.2. Survival data. To increase computational efficiency (particularly true for the scenarios with more individuals captured and a shorter time series), we simulated and fitted the capture-mark-recapture data in the m-array format, using a multinomial likelihood (Burnham, 1987). The data is in the form of two ($T-1 \times T$) matrices $M^J$ and $M^A$, one for each age class generically denoted ($a) \in \{A, J\}$. We have $M^{(a)} = (m_{t,j}^{(a)})$, where $m_{t,j}^{(a)}$ is the number of individuals captured, marked and released in year $t$ that were resighted in year $j + 1$, with $M^{(a)} = (m_{t,j}^{(a)})$, with $m_{t,j}^{(a)} = 0, \forall j < t$. $T$ is the total number of years of capture-recapture history. $m_{t,T}^{(a)}$ is the number of individuals of released at age class ($a$) at time $t$ that were re-sighted the following year, and the last column $m_{t,T}^{(a)}$ is, by convention, the number of individuals released at age class ($a$) at time $t$ that were never re-sighted. We then have:

$$ (4) \qquad \mathbf{m}_{t, \bullet}^{(a)} = (m_{t,t}^{(a)}, m_{t,t+1}^{(a)}, \dots, m_{t,T}^{(a)}) \sim \text{Multinomial} \left( R_t^{(a)}, (\theta_{t,t}^{(a)}, \dots, \theta_{t,T}^{(a)}) \right) $$

with $R_t^{(a)} = \sum_{k=t}^{T} m_{t,k}^{(a)}$ the number of individuals of age class ($a$) released at time $t$.

It is important to note that for the matrix of released juveniles $M^J$, $R_t^J$ corresponds to the number of juveniles newly marked at time $t$. However, $R_t^A$ corresponds to the number of newly marked adults (lets denote it $R_{m,t}^A$), but also of all previously marked juveniles and adults that were released at time $t$. That is,

$$ (5) \qquad R_1^A = R_{m,1}^A, $$
but

$$ (6) \qquad R_{t+1}^A = R_{m,t}^A + \underbrace{\sum_{k=1}^t m_{k,t}^J + m_{k,t}^A}_{\text{previously marked}}. $$

Therefore, unless individuals are not released when marked (e.g., killed, or taken to be released outside of the study population), one needs to provide data and model the number of released adults re-sighted, even if no individuals are first marked as adults. As no individuals are marked as adults here, $R_{m,t}^A = 0$, and so that Equations 5 and 6 can be simplified accordingly.

Note also that in such case were no adults are newly marked, no data on $R_t^A$ is needed to simulate and fit $M^A$. Since $R_t^A = 0$, we have:

$$ (7) \qquad \mathbf{m}_{1,\bullet}^A = 0, $$

$$ (8) \qquad \mathbf{m}_{2,\bullet}^A \sim \text{Multinomial} \left( m_{1,1}^J + m_{1,1}^A, (\theta_{1,t}^A, \dots, \theta_{1,T}^A) \right), $$

$$ (9) \qquad \mathbf{m}_{3,\bullet}^A \sim \text{Multinomial} \left( \sum_{k=1}^2 m_{k,2}^J + m_{k,2}^A, (\theta_{2,t}^A, \dots, \theta_{2,T}^A) \right), $$

and so on.

For juveniles, diagonal elements of the $\theta^J$ matrix write:

$$ \theta_{t,t}^J = \phi_t^J p, $$

with $\phi_t^J$ the first year (i.e. juvenile) survival probability from year $t$ to year $t+1$ (for the species considered), and $p$ the recapture (or re-sighting) probability set as constant among years and age classes, and for $t < j < T$

$$ \theta_{t,j}^J = \phi_t^J \left( \prod_{k=t+1}^j \phi_k^A \right) (1-p)^{j-t} p, $$

with $\phi_t^A$ the adult survival probability from year $t$ to year $t+1$ (for the species considered). $\theta_{t,j}^J$ is the probability of being marked and released as juvenile in year $t$ and recaptured in year $j+1$ as an adult, if $j > t$. The last element pertains to individuals never recaptured

$$ \theta_{t,T}^J = 1 - \sum_{k=t}^{T-1} \theta_{t,k}^J. $$

Similarly for $\theta^A$, the above mentioned equations are identical to the exception that $\phi^J$ is replaced by $\phi^A$, which leads to:

$$ \theta_{t,t}^A = \phi_t^A p, $$

for the diagonal elements of the $\theta^A$ matrix, and for $t < j < T$:

$$ \theta_{t,j}^A = \left( \prod_{k=t}^j \phi_k^A \right) (1-p)^{j-t} p. $$

The last element again pertains to individuals never recaptured

$$ \theta_{t,T}^A = 1 - \sum_{k=t}^{T-1} \theta_{t,k}^A. $$

2.1.3. Fecundity data. Fecundity was modelled using a Poisson regression:

$$ (10) \qquad F_t \sim \text{Poisson}(f_t R_t) $$

with $F_t$ the total number of offspring counted, $R_t$ the number of surveyed broods/litters, and $f_t$ the expected number of offspring (male + female) per adult female each year $t$.
## 2.2. Alternative scenarios and parameter values

### 2.2.1. Density dependence and random temporal variation on demographic rates.
Intra- and inter-species density dependence of survival rates $\phi_{i,t}^{(a)}$ (with $i \in \{V, P\}$) and $(a) \in \{J, A\}$) and fecundities $f_{i,t}$ were modelled on the logit and log scale, respectively. We initially used the same equations as the 2019 model, which are:

$$ (11) \qquad \operatorname{logit}(\phi_{P,t}^J) = \alpha_1 + \alpha_2 n_{P,t}^A $$

where the number of adult predators negatively affects juvenile predator survival (negative intraspecific density dependence),

$$ (12) \qquad \operatorname{logit}(\phi_{V,t}^J) = \alpha_3 + \alpha_4 n_{P,t}^A $$

where the number of adult predators negatively affects juvenile prey survival (predation),

$$ (13) \qquad \operatorname{logit}(\phi_{P,t}^A) = \alpha \phi_P^A $$

$$ (14) \qquad \operatorname{logit}(\phi_{V,t}^A) = \alpha \phi_V^A $$

(no density dependence on adult survival)

$$ (15) \qquad \log(f_{P,t}) = \alpha_5 + \alpha_6 n_{V,t}^J $$

where the number of juvenile prey individuals positively affects predator fecundity, and

$$ (16) \qquad \log(f_{V,t}) = \alpha_7 + \alpha_8 n_{V,t}^A $$

where the number of adult prey individuals negatively affect prey fecundity (negative intraspecific density dependence). Associated results can be found in Supplementary Information Table S1 and Figures S5 to S8.

However, to limit posterior correlation between intercept and slope parameters and improve their estimation, we centered the abundances in the density dependent functions. While centering is typically done and most efficient on mean values, mean abundances varied here from a simulation to the next due to stochasticity. Therefore, intercept parameter values would have to be redefined for each simulation to maintain equivalent mean demographic rate values and asymptotic stage specific abundance equilibria for all simulation. To avoid these complications, we centered by subtracting the corresponding fixed point equilibria estimated in Barraquand and Gimenez (2019) as $\hat{n}_P^A = 21, \hat{n}_V^J = 101$ and $\hat{n}_V^A = 152$. The new $\alpha$ intercept parameters obey the following centered formulas:

$$ (17) \qquad \operatorname{logit}(\phi_{P,t}^J) = \alpha_1 + \alpha_2(n_{P,t}^A - \hat{n}_P^A) $$

$$ (18) \qquad \operatorname{logit}(\phi_{V,t}^J) = \alpha_3 + \alpha_4(n_{P,t}^A - \hat{n}_P^A) $$

$$ (19) \qquad \operatorname{logit}(\phi_{P,t}^A) = \alpha \phi_P^A $$

$$ (20) \qquad \operatorname{logit}(\phi_{V,t}^A) = \alpha \phi_V^A $$

$$ (21) \qquad \log(f_{P,t}) = \alpha_5 + \alpha_6(n_{V,t}^J - \hat{n}_V^J) $$

$$ (22) \qquad \log(f_{V,t}) = \alpha_7 + \alpha_8(n_{V,t}^A - \hat{n}_V^A) $$

To maintain equivalent dynamics to parameter set 1 of the 2019 model, we calculated the intercepts $\alpha_1, \alpha_3, \alpha_5$ and $\alpha_7$ as their original values plus the original slope multiplied by the estimated fixed point equilibrium of the $n$ responsible for density dependence. For example, we now use whenever simulating $\alpha_3 = 0.5 - 0.025 \times 21 = -0.025$ and $\alpha_5 = 0 + 0.004 \times 101 = 0.404$ (Table 1).

In addition, we introduced scenarios with inter-annual random variation in the intercepts of density-dependent links, such that
$$ (23) \qquad \operatorname{logit}(\phi_{P,t}^J) = \alpha_1 + \alpha_2(n_{P,t}^A - \hat{n}_P^A) + \sigma_{\phi_P^J} \epsilon_{\phi_P^J} $$

$$ (24) \qquad \operatorname{logit}(\phi_{V,t}^J) = \alpha_3 + \alpha_4(n_{P,t}^A - \hat{n}_P^A) + \sigma_{\phi_V^J} \epsilon_{\phi_V^J} $$

$$ (25) \qquad \operatorname{logit}(\phi_{P,t}^A) = \alpha_{\phi_P^A} + \sigma_{\phi_P^A} \epsilon_{\phi_P^A} $$

$$ (26) \qquad \operatorname{logit}(\phi_{V,t}^A) = \alpha_{\phi_V^A} + \sigma_{\phi_V^A} \epsilon_{\phi_V^A} $$

with $\epsilon \sim \mathcal{N}(0, 1)$ i.i.d. and

$$ (27) \qquad \log(f_{P,t}) \sim \mathcal{N}(\alpha_5 + \alpha_6(n_{V,t}^J - \hat{n}_V^J), \sigma_{f_P}^2) $$

$$ (28) \qquad \log(f_{V,t}) \sim \mathcal{N}(\alpha_7 + \alpha_8(n_{V,t}^A - \hat{n}_V^A), \sigma_{f_V}^2). $$

Although mathematically identical, we used a parameterisation of the form $\mu + \epsilon\sigma$, $\epsilon \sim \mathcal{N}(0, \sigma^2)$ (sometimes called non-centered) for survival estimates and a centered parameterisation ($\mathcal{N}(\mu, \sigma^2)$) for fecundity estimates as it was found to be optimal for the mixing of the MCMC chains. As we were primarily interested in the ability of multispecies IPMs to estimate species interactions when these were in fact absent, inter species density dependence parameter values for $\alpha_2$ and $\alpha_4$ were either set to zero for the simulations, or at the same value as the 2019 model. Parameter values used to simulate data and their interpretation can be found in Table 1.

<table><caption>**Table 1** – Model parameters with their values. Values of $\alpha_4$ and $\alpha_6$ in the scenarios with true presence of species interactions are presented in parentheses. Temporal standard deviations (SD) are only present in the scenarios with random temporal variation. For interpretation, note that $\alpha_i$ and temporal SD parameters are within exponential functions. For instance, $\alpha_5 = 0.404$ corresponds to a mean fecundity of $e^{0.404} \approx 1.5$.</caption><thead><tr><th>Parameter</th><th>Value</th><th>Interpretation</th></tr></thead><tbody><tr><td>α<sub>1</sub></td><td>0.29</td><td>juvenile predator survival - intercept</td></tr><tr><td>α<sub>2</sub></td><td>-0.01</td><td>juvenile predator survival - slope</td></tr><tr><td>α<sub>3</sub></td><td>-0.025</td><td>juvenile prey survival - intercept</td></tr><tr><td>α<sub>4</sub></td><td><b>0</b> (-0.025)</td><td>juvenile prey survival - slope - inter species density dependence</td></tr><tr><td>α<sub>5</sub></td><td>0.404</td><td>predator fecundity - intercept</td></tr><tr><td>α<sub>6</sub></td><td><b>0</b> (0.004)</td><td>predator fecundity - slope - inter species density dependence</td></tr><tr><td>α<sub>7</sub></td><td>1.24</td><td>prey fecundity - intercept</td></tr><tr><td>α<sub>8</sub></td><td>-0.005</td><td>prey fecundity - slope</td></tr><tr><td>p</td><td>0.7</td><td>recapture probability</td></tr><tr><td>α<sub>φP</sub><sup>A</sup></td><td>logit(0.7)</td><td>adult predator survival - intercept</td></tr><tr><td>α<sub>φV</sub><sup>A</sup></td><td>logit(0.6)</td><td>adult prey survival - intercept</td></tr><tr><td>σ<sub>obs</sub></td><td>0.1</td><td>observation error</td></tr><tr><td>σ<sub>fP</sub></td><td>0.1</td><td>temporal SD of predator fecundity</td></tr><tr><td>σ<sub>fV</sub></td><td>0.1</td><td>temporal SD of prey fecundity</td></tr><tr><td>σ<sub>φP</sub><sup>J</sup></td><td>0.1</td><td>temporal SD of juvenile predator survival</td></tr><tr><td>σ<sub>φP</sub><sup>A</sup></td><td>0.1</td><td>temporal SD deviation of adult predator survival</td></tr><tr><td>σ<sub>φV</sub><sup>J</sup></td><td>0.1</td><td>temporal SD deviation of juvenile prey survival</td></tr><tr><td>σ<sub>φV</sub><sup>A</sup></td><td>0.1</td><td>temporal SD deviation of adult prey survival</td></tr></tbody></table>

2.2.2. Initial values and monitoring setup. For all simulation scenarios in the main text, we used the initial population size vector

$$ \begin{pmatrix} n_{V,1}^J \\ n_{V,1}^A \\ n_{P,1}^J \\ n_{P,1}^A \end{pmatrix} = \begin{pmatrix} 100 \\ 100 \\ 20 \\ 20 \end{pmatrix}, $$
a study period of $T = 30$ years, the yearly number of monitored prey and predator broods/litters respectively $R_t^V = 50$ and $R_t^P = 20$, and the yearly number of marked juveniles was 100 for both species. Results using the monitoring setups of Barraquand and Gimenez (2019) with either 100 marked juveniles per species per year for $T = 10$ years, or 20 marked juveniles per species per year for $T = 30$ years (and the non-centered density-dependencies) are also presented in the Supplementary Information B.

We consider two alternative situations without interspecific interactions: with or without random temporal noise. To compare model performances in the no-interactions setting to cases with interspecific interactions, we also simulated and fitted data in presence of species interactions using the same $\alpha_i$ values as Barraquand and Gimenez (2019) under the four above-mentioned scenarios (i.e., with/without interactions × with/without stochasticity on interactions; see Supplementary Table 2 and Figures S1 and S3 in addition to main text results). For each of these four combinations of parameter sets, we simulated 100 datasets using the Nimble package (de Valpine et al., 2022, 2017, version 0.12.2) in R (R Core Team, 2022, version 4.2.1).

## 2.3. Prior specification and model fitting

Multispecies IPMs were implemented in a Bayesian framework, hence the need to specify priors. When fitting the models to simulated data, we used $N(100, 10)$ and $N(20, 10)$ priors for the initial stage-specific prey and predator population sizes (truncated to be positive). These priors also differed from the 2019 model where they were all set to $N(25, 10^{-5})$.

Priors for standard deviations were chosen as $\sigma \sim \text{Exp}(1)$, which corresponds to priors with maximum entropy on the log and logit scales (e.g., McElreath, 2020). Prior probabilities of recapture were drawn as $p \sim \text{Unif}(0, 1)$ and vital rate/interaction parameters were given weakly informative priors $\alpha_k \sim N(0, 1)$ ($k \in \{1, ..., 8\}$).

Data were both simulated and fitted using the Nimble R package (de Valpine et al., 2022, 2017; R Core Team, 2022, version 0.12.2). To improve their mixing and minimize their posterior correlations, intercepts, slopes and temporal SD were block sampled using automated factor slice samplers (Ponisio et al., 2020; Tibbits et al., 2014). For each simulated dataset, we fitted the same multispecies IPM that was used to generate the data (e.g., no random temporal noise estimated on data without temporal noise), except in that species interactions were estimated even in absence of such interactions. Two MCMC chains were run for 60200 iterations and we sampled the last 60000 iterations every 60th iteration leading to 2000 posterior samples saved per dataset. Real parameter values were used as initial values to minimise time to convergence (see Appendix Section C for an evaluation of the influence of initial values on parameter estimation). We assess convergence and mixing of the chains for all $\alpha_i$ by calculating the potential scale reduction factor ($\hat{R}$, Brooks and Gelman 1998; Gelman and Rubin 1992) and effective sample size ($n_{eff}$) using the "gelman.Diag()" and the "effectiveSize()" functions of the coda package (Plummer et al., 2006, version 0.19-4). We only used outputs from models for which all $\alpha_i$ had $\hat{R} < 1.1$ and $n_{eff} > 50$, that is, 100/100 models for the scenario without random temporal variation and 94/100 models for the scenario with random temporal variation. The computer code is provided at https://github.com/MatthieuPaquet/multi物种.

## 3. Results

Overall, estimates of density dependence curves were unbiased, regarding interspecific density dependence (either absent, Figures 1 and 3, or present, Figures S1 and S3) as well as intraspecific density dependence. This was true without and with temporal stochasticity (Figures 1 to 4).

This absence of bias extends to the alternative data designs with smaller sample sizes considered in Barraquand and Gimenez (2019) (shown in Supplementary Information in Figures S5 to S8). Estimated $\alpha_i$ parameters also did not show sign of bias in any scenario (Table 2 and Table S1).
<figure><img src="image_3.png" /></figure>

**Figure 1** – Density-dependencies for juvenile survival rates (A for predator and B for prey) as well as prey (C) and predator (D) fecundities in the scenario without random time variation. Purple: simulated relationships, light green: posterior mean relationships for all 100 fitted models, dark green: average of the posterior mean relationships. True inter species density-dependencies (right panels) were set to be absent.
<figure><img src="image_4.png" /></figure>

**Figure 2** – Example of posterior mean (blue-green line) and 95% Credible Intervals (grey polygons) of density-dependencies for juvenile survival rates (A for predator and B for prey) as well as prey (C) and predator (D) fecundities estimated by one of the 100 models run in the scenario without random time variation. Purple lines indicate the simulated (true) relationships. Points represent estimated mean demographic parameter each year plotted against estimated yearly abundance values, and vertical and horizontal error bars their respective 95% Credible Intervals.
<figure><img src="image_5.png" /></figure>

**Figure 3** – Density-dependencies for juvenile survival rates (A for predator and B for prey) as well as prey (C) and predator (D) fecundities in the scenario with random time variation. Purple: simulated relationships, light green: posterior mean relationships for the 94 fitted models that appear to converge satisfactorily, dark green: average of the posterior mean relationships. True inter species density-dependencies (right panels) were set to be absent.
<figure><img src="image_6.png" /><figcaption>**Figure 4** – Example of posterior mean (blue-green line) and 95% Credible Intervals (grey polygons) of density-dependencies for juvenile survival rates (A for predator and B for prey) as well as prey (C) and predator (D) fecundities estimated by one of the 100 models run in the scenario with random time variation. Purple lines indicate the simulated (true) relationships. Points represent estimated mean demographic parameter each year plotted against estimated yearly abundance values, and vertical and horizontal error bars their respective 95% Credible Intervals.</figcaption></figure>

We did not detect more false positive species interactions than expected by chance when investigating the coverage probability of the species interaction parameters at 95% (i.e., the proportion of simulations where 95% Crl of estimated parameter includes the true parameter value). In the scenario with 100 juveniles marked each year for 30 years and no interspecific density dependence nor temporal random variation, this probability was 0.95 for $\alpha_4$ and 0.92 for $\alpha_6$ (cf Table 2, see Figure 2 for an example of estimated mean and pointwise 95% Crl density dependent curves). Coverage probabilities were also satisfactory when interspecific interactions were simulated to be nonzero (0.94 and 0.93). Species interactions parameters were still estimable with no noticeable bias in the presence of random time variation (Figures 3 and 4), in which case the coverage probabilities of the species interaction parameters $\alpha_4$ and $\alpha_6$ at 95% were 0.99 and 0.98 respectively in absence of interspecific interactions (Table 2). In the presence of interspecific interactions, coverage values were both 0.96. Moreover, the addition of random
time variation did not noticeably alter the precision of the species interaction parameters, both in absence and presence of species interactions (Figure S3, Table 2).

<table><caption>**Table 2** – Summary table of parameter estimates. Value refers to the true values used to simulate the data and values of the interspecific density dependent parameters are highlighted in bold. Estimate (95% quantiles) are the mean and the 95% quantiles of the posterior mean estimates. Coverage 95% is the proportion of 95% Credible Intervals that included the true parameter values.</caption><thead><tr><th>Scenario</th><th>Param.</th><th>Value</th><th>Estimate (95% quantiles)</th><th>Coverage 95%</th></tr></thead><tbody><tr><td>30 years</td><td>α<sub>1</sub></td><td>0.29</td><td>0.304 (0.166; 0.482)</td><td>0.97</td></tr><tr><td>100 ind. marked/year</td><td>α<sub>2</sub></td><td>-0.01</td><td>-0.013 (-0.032; 0.005)</td><td>0.95</td></tr><tr><td>No temporal noise</td><td>α<sub>3</sub></td><td>-0.025</td><td>-0.033 (-0.144; 0.083)</td><td>0.99</td></tr><tr><td><strong>No interspecies DD</strong></td><td>α<sub>4</sub></td><td><strong>0</strong></td><td>0.001 (-0.017; 0.019)</td><td>0.95</td></tr><tr><td></td><td>α<sub>5</sub></td><td>0.404</td><td>0.413 (0.199; 0.635)</td><td>0.93</td></tr><tr><td></td><td>α<sub>6</sub></td><td><strong>0</strong></td><td>0 (-0.008; 0.007)</td><td>0.92</td></tr><tr><td></td><td>α<sub>7</sub></td><td>1.24</td><td>1.243 (1.198; 1.287)</td><td>0.97</td></tr><tr><td></td><td>α<sub>8</sub></td><td>-0.005</td><td>-0.005 (-0.006; -0.004)</td><td>0.96</td></tr><tr><td>30 years</td><td>α<sub>1</sub></td><td>0.29</td><td>0.281 (0.124; 0.428)</td><td>0.947</td></tr><tr><td>100 ind. marked/year</td><td>α<sub>2</sub></td><td>-0.01</td><td>-0.013 (-0.038; 0.003)</td><td>0.947</td></tr><tr><td>Temporal noise</td><td>α<sub>3</sub></td><td>-0.025</td><td>-0.02 (-0.173; 0.122)</td><td>0.947</td></tr><tr><td><strong>No interspecies DD</strong></td><td>α<sub>4</sub></td><td><strong>0</strong></td><td>0.001 (-0.012; 0.017)</td><td>0.989</td></tr><tr><td></td><td>α<sub>5</sub></td><td>0.404</td><td>0.399 (0.236; 0.532)</td><td>0.968</td></tr><tr><td></td><td>α<sub>6</sub></td><td><strong>0</strong></td><td>0 (-0.005; 0.004)</td><td>0.979</td></tr><tr><td></td><td>α<sub>7</sub></td><td>1.24</td><td>1.244 (1.17; 1.319)</td><td>0.968</td></tr><tr><td></td><td>α<sub>8</sub></td><td>-0.005</td><td>-0.005 (-0.007; -0.004)</td><td>0.968</td></tr><tr><td>30 years</td><td>α<sub>1</sub></td><td>0.29</td><td>0.282 (0.122; 0.427)</td><td>0.97</td></tr><tr><td>100 ind. marked/year</td><td>α<sub>2</sub></td><td>-0.01</td><td>-0.011 (-0.027; 0.004)</td><td>0.98</td></tr><tr><td>No temporal noise</td><td>α<sub>3</sub></td><td>-0.025</td><td>-0.019 (-0.156; 0.125)</td><td>0.99</td></tr><tr><td><strong>Interspecies DD</strong></td><td>α<sub>4</sub></td><td><strong>-0.025</strong></td><td>-0.026 (-0.042; -0.009)</td><td>0.94</td></tr><tr><td></td><td>α<sub>5</sub></td><td>0.404</td><td>0.395 (0.264; 0.502)</td><td>0.94</td></tr><tr><td></td><td>α<sub>6</sub></td><td><strong>0.004</strong></td><td>0.004 (-0.001; 0.01)</td><td>0.93</td></tr><tr><td></td><td>α<sub>7</sub></td><td>1.24</td><td>1.241 (1.195; 1.281)</td><td>0.95</td></tr><tr><td></td><td>α<sub>8</sub></td><td>-0.005</td><td>-0.005 (-0.006; -0.004)</td><td>0.96</td></tr><tr><td>30 years</td><td>α<sub>1</sub></td><td>0.29</td><td>0.29 (0.143; 0.459)</td><td>0.967</td></tr><tr><td>100 ind. marked/year</td><td>α<sub>2</sub></td><td>-0.01</td><td>-0.01 (-0.026; 0.007)</td><td>0.989</td></tr><tr><td>Temporal noise</td><td>α<sub>3</sub></td><td>-0.025</td><td>-0.014 (-0.176; 0.145)</td><td>0.967</td></tr><tr><td><strong>Interspecies DD</strong></td><td>α<sub>4</sub></td><td><strong>-0.025</strong></td><td>-0.025 (-0.041; -0.007)</td><td>0.957</td></tr><tr><td></td><td>α<sub>5</sub></td><td>0.404</td><td>0.403 (0.276; 0.523)</td><td>0.957</td></tr><tr><td></td><td>α<sub>6</sub></td><td><strong>0.004</strong></td><td>0.004 (-0.001; 0.009)</td><td>0.957</td></tr><tr><td></td><td>α<sub>7</sub></td><td>1.24</td><td>1.237 (1.179; 1.311)</td><td>0.913</td></tr><tr><td></td><td>α<sub>8</sub></td><td>-0.005</td><td>-0.005 (-0.008; -0.004)</td><td>0.924</td></tr></tbody></table>

## 4. Discussion

Building on the multispecies integrated predator-prey model of Barraquand and Gimenez (2019), we investigated here whether multispecies IPMs could output interactions where there are in fact none. We did so by modelling functions relating vital rates to stage-specific species densities, whose slope parameters are used to model species interactions. We found that when those slopes were simulated as zero, the estimates were centered on zero and therefore unbiased. There was also a good coverage probability of interaction parameters (close to 0.95 for 95% CrIs). We also found that adding temporal variability to these multispecies density-dependent link functions did not alter these results. This confirms that multispecies IPMs are a promising way to estimate species interactions, and in particular, that they could be used to infer whether two species interact or not when such information is missing.
These results are encouraging, though some readers might find our sample sizes relatively large (see Appendix B for slightly lower sample sizes). In a previous version of this work, we inadvertently omitted the $\theta^A$ CMR array in the code, which transformed the model into a capture-removal model (i.e., individuals were re-captured only once and then removed from the population, as in hunting or fishing data). In this configuration, the lower amount of data on survival and detection provided proper estimation of all quantities for the main text data design but not those of Appendix B for which convergence was not always reached. With live capture-recapture data, all data designs (main text and Appendix B) now provide satisfactory convergence and estimation. Moreover, in field population studies, additional types of data available are likely to improve the estimation of species interactions and we give three examples below. First, when age classes can be determined during the count observation process, including such information explicitly in the model (see e.g., Paquet et al., 2019; Weegman et al., 2016) will increase identifiability and/or precision of survival parameters and age specific abundances, and therefore will likely improve the estimation of density dependence parameters as well. This stage-specific abundance information may also allow, in some cases where counts are provided with little error, to remove the observation process, which we cannot do in our current model formulation because the observed population size sums adult and juvenile densities, and this sum has to arise from a probability distribution (Equations (2) and (3)). Second, integrating dead prey recovery data is likely to give extra information on the strength of predator-prey interactions. Dead recoveries are classically implemented in capture-mark-recovery models (North and Morgan, 1979; Seber, 1972) which in some cases can be combined with CMR data (Barker, 1999) and counts (Reynolds et al., 2009). Since the probability to find a dead prey is likely affected by predation rates in the population (e.g., in some systems prey eaten will not be recovered, in others dead recoveries may present signs of predation), taking the predation process into account in the dead recoveries data-generation mechanism could improve the estimation of the strength of predator-prey interactions. Finally, the spatial structure of the data should contain additional information that may help to estimate parameters. The extension to spatially explicit IPMs (Chandler and Clark, 2014; Zhao, 2020) for interacting populations represents a promising way forward for the estimation of species interactions.

We commented above on the amount of data and possible additional data types. However, the efficiency of multispecies IPMs in estimating species interactions may also depend on the parameter set, and thus on the ecological features of the populations studied. For example, the parameters considered here correspond well to vertebrate predator-prey systems with a stable equilibrium in absence of environmental perturbations. Faster life histories, different stage or age structure, and multiple factors contributing to altering the quantity of information encapsulated in the various data streams may alter the sample sizes required for efficient inferences. When applying these models to new systems with different life history parameters and density-dependent structures (e.g., predators also eating adult prey), simulated datasets with plausible ecological features for the empirical system considered (and similar data designs), will help confirm that parameter values can be recovered without bias and with sufficient precision. Tools such as JAGS (Plummer, 2003) or Nimble (de Valpine et al., 2017) make it particularly handy to both simulate and fit data with complex dynamic models.

Finally, while using the same model to simulate and fit the data is a necessary first step to (i) assess the identifiability of model parameters (and assess the amount and type of data needed for practical identifiability), (ii) evaluate the coverage of parameter estimates, and (iii) check for bias in the estimates that can still occur, notably because of limited sample sizes (Paquet et al., 2021), an important next step will be to evaluate the sensitivity of multi-species IPM estimates to model mis-specifications (Plard et al., 2021). For example, different functions than the log and logit links chosen here may be used to fit or to simulate intra- and inter-specific density-dependencies. Hence, we encourage future work to try fitting a broader range of plausible models that differ from the model used to simulate the data (or conversely, to simulate from more mechanistic models) in order to assess such sensitivity.
## Acknowledgements

We thank Olivier Gimenez for discussion of these analyses and comments on the manuscript, as well as the reviewers and recommenders for their constructive feedback. Colleagues from IMB and Biogeco labs are thanked for a welcoming environment. Preprint version 2 of this article has been peer-reviewed and recommended by Peer Community In Ecology (https://doi.org/10.24072/pci.ecology.100522, Coulson and Alonso (2023)).

## Data and code

Computer R code used to generate and analyze the data is available at https://github.com/MatthieuPaquet/multi物种 and at https://doi.org/10.17605/OSF.IT/XFA6E (Paquet and Barraquand, 2023) together with the generated data and figures.

## Funding

Funding was provided through grant ANR-20-CE45-0004 to FB.

## Conflict of interest disclosure

The authors of this preprint declare that they have no financial conflict of interest with respect to the content of this article. Matthieu Paquet and Frédéric Barraquand are recommenders at PCI Ecology.

## References

* Abadi F, Gimenez O, Ullrich B, Arlettaz R, Schaub M (2010). Estimation of immigration rate using integrated population models. *Journal of Applied Ecology* **47**, 393-400. https://doi.org/10.1111/j.1365-2664.2010.01789.x.

* Auger-Méthé M, Field C, Albertsen CM, Derocher AE, Lewis MA, Jonsen ID, Mills Flemming J (2016). State-space models' dirty little secrets: even simple linear Gaussian models can have estimation problems. *Scientific reports* **6**, 1-10. https://doi.org/10.1038/srep26677.

* Barker RJ (1999). *Joint analysis of mark-recapture, resighting and ring-recovery data with age-dependence and marking-effect*. *Bird Study* **46**, S82-S91. https://doi.org/10.1080/00063659909477235.

* Barraquand F, Gimenez O (2019). *Integrating multiple data sources to fit matrix population models for interacting species*. *Ecological modelling* **411**, 1087-113. https://doi.org/10.1016/j.ecolmodel.2019.06.001.

* Barraquand F, Picoche C, Detto M, Hartig F (2021). *Inferring species interactions using Granger causality and convergent cross mapping*. *Theoretical Ecology* **14**, 87-105. https://doi.org/10.1007/s12080-020-00482-7.

* Besbeas P, Freeman SN, Morgan BJT, Catchpole EA (2002). *Integrating mark-recapture-recovery and census data to estimate animal abundance and demographic parameters*. *Biometrics* **58**, 540-547. https://doi.org/10.1111/j.0006-341X.2002.00540.x.

* Brooks SP, Gelman A (1998). *General methods for monitoring convergence of iterative simulations*. *Journal of computational and graphical statistics* **7**, 434-455. https://doi.org/10.1080/10618600.1998.10474787.

* Burnham KP (1987). *Design and analysis methods for fish survival experiments based on release-recapture*. American Fisheries Society.

* Chandler RB, Clark JD (2014). *Spatially explicit integrated population models*. *Methods in Ecology and Evolution* **5**, 1351-1360. https://doi.org/10.1111/2041-210X.12153.

* Coulson T, Alonso D (2023). *Addressing the daunting challenge of estimating species interactions from count data*. *Peer Community in Ecology* **1**, 100522. https://doi.org/10.24072/pci.ecology.100522.
de Valpine P, Paciorek C, Turek D, Michaud N, Anderson-Bergman C, Obermeyer F, Wehrhahn Cortes C, Rodríguez A, Temple Lang D, Paganin S, Hug J (2022). NIMBLE: MCMC, Particle Filtering, and Programmable Hierarchical Modeling. Version 0.12.2. https://doi.org/10.5281/zenodo.1211190. URL: https://cran.r-project.org/package=nimble.

de Valpine P, Turek D, Paciorek C, Anderson-Bergman C, Temple Lang D, Bodik R (2017). Programming with models: writing statistical algorithms for general model structures with NIMBLE. Journal of Computational and Graphical Statistics 26, 403-417. https://doi.org/10.1080/10618600.2016.1172487.

Dennis B, Ponciano JM, Lele SR, Taper ML, Staples DF (2006). Estimating density dependence, process noise, and observation error. Ecological Monographs 76, 323-341. https://doi.org/10.1890/0012-9615(2006)76[323:EDDPNA]2.0.CO;2.

Gelman A, Rubin DB (1992). Inference from iterative simulation using multiple sequences. Statistical science, 457-472. https://doi.org/10.1214/ss/1177011136.

Kéry M, Schaub M (2011). Bayesian population analysis using WinBUGS: a hierarchical perspective. Academic Press. https://doi.org/10.1016/b978-0-12-387020-9.00001-8.

Knape J (2008). Estimability of density dependence in models of time series data. Ecology 89, 2994-3000. https://doi.org/10.1890/08-0071.1.

Lahoz-Monfort JJ, Harris MP, Wanless S, Freeman SN, Morgan BJT (2017). Bringing it all together: multi-species integrated population modelling of a breeding community. Journal of Agricultural, Biological and Environmental Statistics 22, 140-160. https://doi.org/10.1007/s13253-017-0279-4.

McElreath R (2020). Statistical rethinking: A Bayesian course with examples in R and Stan. Chapman and Hall/CRC. https://doi.org/10.1201/9780429029608.

Mutshinda CM, O'Hara RB, Woiwod IP (2009). What drives community dynamics? Proceedings of the Royal Society B: Biological Sciences 276, 2923-2929. https://doi.org/10.1098/rspb.2009.0523.

North PM, Morgan BJT (1979). Modelling heron survival using weather data. Biometrics, 667-681. https://doi.org/10.2307/2530260.

Paquet M, Arlt D, Knape J, Low M, Forslund P, Pärt T (2019). Quantifying the links between land use and population growth rate in a declining farmland bird. Ecology and evolution 9, 868-879. https://doi.org/10.1002/ece3.4766.

Paquet M, Barraquand F (2023). Data and code: Assessing species interactions using integrated predator-prey models. https://doi.org/10.17605/OSF.IO/XFA6E.

Paquet M, Knape J, Arlt D, Forslund P, Pärt T, Flagstad Ø, Jones CG, Nicoll MAC, Norris K, Pemberton JM, Sand H, Svensson L, Tatayah V, Wabakken P, Wikenros C, Åkesson M, Low M (2021). Integrated population models poorly estimate the demographic contribution of immigration. Methods in Ecology and Evolution 12, 1899-1910. https://doi.org/10.1111/2041-210X.13667.

Péron G, Koons DN (2012). Integrated modeling of communities: parasitism, competition, and demographic synchrony in sympatric ducks. Ecology 93, 2456-2464. https://doi.org/10.1890/11-1881.1.

Plard F, Turek D, Schaub M (2021). Consequences of violating assumptions of integrated population models on parameter estimates. Environmental and Ecological Statistics 28, 667-695. https://doi.org/10.1007/s10651-021-00507-2.

Plummer M (2003). JAGS: A program for analysis of Bayesian graphical models using Gibbs sampling. In: Proceedings of the 3rd international workshop on distributed statistical computing. Vol. 124. 125.10. Vienna, Austria., pp. 1-10.

Plummer M, Best N, Cowles K, Vines K (2006). CODA: Convergence Diagnosis and Output Analysis for MCMC. R News 6, 7-11. URL: https://journal.r-project.org/articles/RN-2006-002/RN-2006-002.pdf.

Ponisio LC, Valpine P, Michaud N, Turek D (2020). One size does not fit all: Customizing MCMC methods for hierarchical models using NIMBLE. Ecology and evolution 10, 2385-2416. https://doi.org/10.1002/ece3.6053.
Quéroué M, Barbraud C, Barraquand F, Turek D, Delord K, Pacoureau N, Gimenez O (2021). *Multispecies integrated population model reveals bottom-up dynamics in a seabird predator-prey system*. *Ecological monographs* **91**, e01459. https://doi.org/10.1002/ecm.1459

R Core Team (2022). *R: A Language and Environment for Statistical Computing*. R Foundation for Statistical Computing. Vienna, Austria. URL: https://www.R-project.org/.

Reynolds TJ, King R, Harwood J, Frederiksen M, Harris MP, Wanless S (2009). *Integrated data analysis in the presence of emigration and mark loss*. *Journal of Agricultural, Biological, and Environmental Statistics* **14**, 411–431. https://doi.org/10.1198/jabes.2009.08008.

Sander EL, Wootton JT, Allesina S (2017). *Ecological network inference from long-term presence-absence data*. *Scientific reports* **7**, 1–12. https://doi.org/10.1038/s41598-017-07009-x.

Seber G (1972). *Estimating survival rates from bird-band returns*. *The Journal of Wildlife Management*, 405-413. https://doi.org/10.2307/3799070.

Strydom T, Catchen MD, Banville F, Caron D, Dansereau G, Desjardins-Proulx P, Forero-Muñoz NR, Higino G, Mercier B, Gonzalez A, Gravel D, Pollock L, Poisot T (2021). *A roadmap towards predicting species interaction networks (across space and time)*. *Philosophical Transactions of the Royal Society B* **376**, 20210063. https://doi.org/10.1098/rstb.2021.0063.

Tibbits MM, Groendyke C, Haran M, Liechty JC (2014). *Automated factor slice sampling*. *Journal of Computational and Graphical Statistics* **23**, 543–563. https://doi.org/10.1080/10618600.2013.791193.

Weegman MD, Bearhop S, Fox AD, Hilton GM, Walsh AJ, McDonald JL, Hodgson DJ (2016). *Integrated population modelling reveals a perceived source to be a cryptic sink*. *Journal of Animal Ecology* **85**, 467–475. https://doi.org/10.1111/1365-2656.12481.

Zhao Q (2020). *On the sampling design of spatially explicit integrated population models*. *Methods in Ecology and Evolution* **11**, 1207–1220. https://doi.org/10.1111/2041-210X.13457.
**Supplementary Information**

A. Results for the scenarios with species interactions

<figure><img src="image_7.png" /></figure>

**Figure S1** – Density-dependencies for juvenile survival rates (A for predator and B for prey) as well as prey (C) and predator (D) fecundities in the scenario without random time variation in presence of true inter species density-dependencies. Purple: simulated relationships, light green: posterior mean relationships for the 100 fitted models that appear to converge satisfactorily, dark green: average of the posterior mean relationships.
<figure><img src="image_8.png" /></figure>

**Figure S2** – Example of posterior mean (blue-green line) and 95% Credible Intervals (grey polygons) of density-dependencies for juvenile survival rates (A for predator and B for prey) as well as prey (C) and predator (D) fecundities estimated by one of the 100 models run in the scenario without random time variation in presence of true inter species density-dependencies. Purple lines indicate the simulated (true) relationships. Points represent estimated mean demographic parameter each year plotted against estimated yearly abundance values, and vertical and horizontal error bars their respective 95% Credible Intervals.
<figure><img src="image_9.png" /></figure>

**Figure S3** – Density-dependencies for juvenile survival rates (A for predator and B for prey) as well as prey (C) and predator (D) fecundities in the scenario with random time variation in presence of true inter species density-dependencies. Purple: simulated relationships, light green: posterior mean relationships for the 92 fitted models that appear to converge satisfactorily, dark green: average of the posterior mean relationships.
<figure><img src="image_10.png" /></figure>

**Figure S4** – Example of posterior mean (blue-green line) and 95% Credible Intervals (grey polygons) of density-dependencies for juvenile survival rates (A for predator and B for prey) as well as prey (C) and predator (D) fecundities estimated by one of the 100 models run in the scenario with random time variation in presence of true inter species density-dependencies. Purple lines indicate the simulated (true) relationships. Points represent estimated mean demographic parameter each year plotted against estimated yearly abundance values, and vertical and horizontal error bars their respective 95% Credible Intervals.
## B. Results for the scenarios with 100 juveniles of each species marked each year for 10 years, and 20 juveniles of each species marked for 30 years, without centering abundances in the link functions

The results presented below follow the data design of Barraquand and Gimenez (2019).

<table><caption>**Table S1** – Value refers to the true values used to simulate the data and values of the interspecific density dependent parameters are highlighted in bold. Estimate (95% quantiles) are the mean and the 95% quantiles of the posterior mean estimates. Coverage 95% is the proportion of 95% Credible Intervals that included the true parameter values.</caption><thead><tr><th>Scenario</th><th>Param.</th><th>Value</th><th>Estimate (95% quantiles)</th><th>Coverage 95%</th></tr></thead><tbody><tr><td>10 years</td><td>α<sub>1</sub></td><td>0.5</td><td>0.458 (-0.545; 1.293)</td><td>0.98</td></tr><tr><td>100 ind. marked/year</td><td>α<sub>2</sub></td><td>-0.01</td><td>-0.008 (-0.041; 0.031)</td><td>0.99</td></tr><tr><td rowspan="7">No temporal noise</td><td>α<sub>3</sub></td><td>-0.025</td><td>-0.045 (-0.78; 0.563)</td><td>0.99</td></tr><tr><td>α<sub>4</sub></td><td>0</td><td>0.001 (-0.026; 0.033)</td><td>0.98</td></tr><tr><td>α<sub>5</sub></td><td>0.404</td><td>0.286 (-0.577; 1.098)</td><td>0.949</td></tr><tr><td>α<sub>6</sub></td><td>0</td><td>0.001 (-0.005; 0.008)</td><td>0.96</td></tr><tr><td>α<sub>7</sub></td><td>2</td><td>1.998 (1.773; 2.227)</td><td>0.98</td></tr><tr><td>α<sub>8</sub></td><td>-0.005</td><td>-0.005 (-0.006; -0.004)</td><td>0.97</td></tr><tr><td>α<sub>1</sub></td><td>0.5</td><td>0.278 (-0.361; 0.874)</td><td>0.99</td></tr><tr><td>100 ind. marked/year</td><td>α<sub>2</sub></td><td>-0.01</td><td>-0.001 (-0.028; 0.029)</td><td>0.99</td></tr><tr><td rowspan="7">Temporal noise</td><td>α<sub>3</sub></td><td>-0.025</td><td>-0.023 (-0.594; 0.614)</td><td>1</td></tr><tr><td>α<sub>4</sub></td><td>0</td><td>0 (-0.025; 0.023)</td><td>1</td></tr><tr><td>α<sub>5</sub></td><td>0.404</td><td>0.323 (-0.379; 0.918)</td><td>0.99</td></tr><tr><td>α<sub>6</sub></td><td>0</td><td>0.001 (-0.004; 0.006)</td><td>1</td></tr><tr><td>α<sub>7</sub></td><td>2</td><td>1.911 (1.467; 2.323)</td><td>0.948</td></tr><tr><td>α<sub>8</sub></td><td>-0.005</td><td>-0.004 (-0.007; -0.002)</td><td>0.927</td></tr><tr><td>30 years</td><td>α<sub>1</sub></td><td>0.5</td><td>0.56 (-0.155; 1.278)</td><td>0.99</td></tr><tr><td>20 ind. marked/year</td><td>α<sub>2</sub></td><td>-0.01</td><td>-0.013 (-0.045; 0.021)</td><td>0.99</td></tr><tr><td rowspan="7">No temporal noise</td><td>α<sub>3</sub></td><td>-0.025</td><td>-0.031 (-0.529; 0.436)</td><td>0.97</td></tr><tr><td>α<sub>4</sub></td><td>0</td><td>0.001 (-0.024; 0.023)</td><td>0.97</td></tr><tr><td>α<sub>5</sub></td><td>0.404</td><td>0.329 (-0.305; 1.009)</td><td>0.97</td></tr><tr><td>α<sub>6</sub></td><td>0</td><td>0.001 (-0.005; 0.006)</td><td>0.98</td></tr><tr><td>α<sub>7</sub></td><td>2</td><td>2.009 (1.832; 2.194)</td><td>0.96</td></tr><tr><td>α<sub>8</sub></td><td>-0.005</td><td>-0.005 (-0.006; -0.004)</td><td>0.97</td></tr><tr><td>α<sub>1</sub></td><td>0.5</td><td>0.527 (-0.275; 1.206)</td><td>0.968</td></tr><tr><td>30 years</td><td>α<sub>2</sub></td><td>-0.01</td><td>-0.011 (-0.041; 0.019)</td><td>0.968</td></tr><tr><td rowspan="7">Temporal noise</td><td>α<sub>3</sub></td><td>-0.025</td><td>-0.007 (-0.528; 0.533)</td><td>0.979</td></tr><tr><td>α<sub>4</sub></td><td>0</td><td>0 (-0.02; 0.017)</td><td>0.989</td></tr><tr><td>α<sub>5</sub></td><td>0.404</td><td>0.316 (-0.191; 0.823)</td><td>0.968</td></tr><tr><td>α<sub>6</sub></td><td>0</td><td>0.001 (-0.004; 0.004)</td><td>0.968</td></tr><tr><td>α<sub>7</sub></td><td>2</td><td>1.954 (1.702; 2.247)</td><td>0.957</td></tr><tr><td>α<sub>8</sub></td><td>-0.005</td><td>-0.005 (-0.006; -0.003)</td><td>0.957</td></tr></tbody></table>
<figure><img src="image_11.png" /></figure>

**Figure S5** – Density-dependencies for juvenile survival rates (A for predator and B for prey) as well as prey (C) and predator (D) fecundities in the scenario with 100 juveniles per species marked each year for 10 years without random time variation in absence of true inter species density-dependencies. Purple: simulated relationships, light green: posterior mean relationships for the 99 fitted models that appear to converge satisfactorily, dark green: average of the posterior mean relationships.
<figure><img src="image_12.png" /></figure>

**Figure S6** – Density-dependencies for juvenile survival rates (A for predator and B for prey) as well as prey (C) and predator (D) fecundities in the scenario with 100 juveniles per species marked each year for 10 years with random time variation in absence of true inter species density-dependencies. Purple: simulated relationships, light green: posterior mean relationships for the 96 fitted models that appear to converge satisfactorily, dark green: average of the posterior mean relationships.
<figure><img src="image_13.png" /></figure>

**Figure S7** – Density-dependencies for juvenile survival rates (A for predator and B for prey) as well as prey (C) and predator (D) fecundities in the scenario with 20 juveniles per species marked each year for 30 years without random time variation in absence of true inter species density-dependencies. Purple: simulated relationships, light green: posterior mean relationships for the 100 fitted models that appear to converge satisfactorily, dark green: average of the posterior mean relationships.
<figure><img src="image_14.png" /></figure>

**Figure S8** – Density-dependencies for juvenile survival rates (A for predator and B for prey) as well as prey (C) and predator (D) fecundities in the scenario with 20 juveniles per species marked each year for 30 years with random time variation in absence of true inter species density-dependencies. Purple: simulated relationships, light green: posterior mean relationships for the 94 fitted models that appear to converge satisfactorily, dark green: average of the posterior mean relationships.
### C. Sensitivity of parameter estimation to the choice of initial values

To assess whether the accuracy of the estimation of density dependent parameters was con-
ditioned by the fact that we used true parameter values as initial values, we also ran the MCMC
using values that substantially deviated from the true value and expected posterior distributions.
For this study, we used data (and the corresponding model) without temporal random noise and
without true interspecific interactions. We chose one simulated dataset for which the true val-
ues of α₂, α₄, α₆ and α₈ fell well within the 95% credible intervals of the posterior samples when
using the true value as initial value (see script https://github.com/MatthieuPaquet/multi-
species/blob/main/script_initial_values.R for more details on the procedure). We then
simulated 100 sets of initial values that deviated from the true values by 4 standard deviations
estimated from the posterior samples when the true values were used as initial values (hereafter
$SD_{\hat{\alpha}_i}$). For the parameters for which negative density dependence was expected, we simulated
the 100 initial values as $\alpha_i^{init} \sim N(\alpha_i - 4SD_{\hat{\alpha}_i}, SD_{\hat{\alpha}_i})$ whereas for $\alpha_8$, which was a potentially
positive prey $\rightarrow$ predator link (and would have been assumed positive in an empirical analysis),
we used $\alpha_8^{init} \sim N(\alpha_8 + 4SD_{\hat{\alpha}_8}, SD_{\hat{\alpha}_8})$. We used true parameter values as initial values for all
other model parameters. Preliminary runs showed that convergence was reached very quickly
(typically after a couple of iterations) with efficient mixing. We then ran 2 chains for 1200 it-
erations and discarded the first 200 as burn-in and did not use thinning. For comparison we
also run 2 MCMC chains once, under the same settings, using the true values as initial values
(see script https://github.com/MatthieuPaquet/multi species/blob/main/script_MCMC_\n simulatedinitial\_values\_out\_of\_posterior.R). The results showed no sign of influence of
the initial value chosen on the parameter estimates (Figure S9).
<figure><img src="image_15.png" /></figure>

**Figure S9** – Estimation of density dependent parameter values ($\alpha_2$ in panel A, $\alpha_4$ in panel B, $\alpha_6$ in panel C and $\alpha_8$ in panel D) in relation to the initial values chosen to start the MCMC chains. Dots show the posterior means and vertical lines the 95% credible intervals. Purple horizontal lines highlight the value used to simulate the data. Red dots and intervals show the case where the true values are used as initial values.