Beyond pairwise species interactions: coarser inference of their joined effects is more relevant

This manuscript explores the limits of accurately inferring biotic interaction strength and direction from datasets that are simulated to mimic certain field conditions. The authors answered a common logistical question in community ecology: when can we infer more from sparse data? They provided clear results showing that inference will be more challenging when species are not dispersal limited, but even then the mean and/or variance of the inferred interaction matrix may still be reliable.

I have found the paper very clearly written, and the figures are well-annotated and relatively accessible to a general audience. The Methods text may need a bit more clarity to guide us through the steps (e.g., at times I was not sure how inference was made / how models were fitted on the simulated data; do the mathematical equations describe either or both the simulation and inference?), but there was nothing major that prevents a reader like me from grasping the full picture.

Here are some specific comments:

For a reader like me who is a frequent user of joint species distribution model (or think in that framework), several of the cited papers in this manuscript keep leaving me a lingering thought "so what does this mean to my JSDM or co-distribution datasets?" The authors did quite a good job outlining JSDMs in paragraph 3, as well as their shortcomings citing Poggiato et al. 2021, but in the Discussion I find myself still hungry for more insights from the authors. From here I learned that dispersal limitation (or the lack thereof) is key to whether the results from JSDMs are reliable in inferring biotic interaction, but which results is it from a JSDM? As a user of JSDMs, it is still unclear to me whether we could reliable interpret the latent associations (or residual components) as biotic interactions / associations not accounted for by the environment (assuming that we have measured all the important abiotic covariates), or do we include neighbour densities as measured covariates. Knowing Leibold et al. (2022) Oikos, I was waiting for a bit more insights into the way forward to use JSDMs. Of course, not everything needs to be analysed as a JSDM (or using existing JSDM softwares). Maybe this is why the authors opted not be too overprescriptive. It is just that the regular JSDM citations leave one hanging about them. I apologise for not having any concrete solution to suggest at this moment.

The second point is about the term non-additivity sprikled throughtout this work. I appreciate the authors noting this, but it isn't clear to me whether they think non-addivity is something to avoid, or something to tackle. In the 3rd Discussion paragraph, I learned that non-addivities (or indirect effect) may be more prevalent in the environment-tracking scenario or when most species have reached equilibrium, so do well avoid ill-posed questions under these situations? Or would the inclusion of N_j by E(x) interaction term in a statistical model be the first step to account for non-additivities? I am not expecting a lengthy response but only exploring if the authors have clearer suggestions for this, as I understand that this may not fit within the scope of the current work.

I appreciate the effort going into explaining complex terminologies and concepts, but there are still a few places that could be more spelled out to help us understand. For example:

- In the 5th Intro paragraph, you may need to guide us through direct to net interactions. Why does the inverse tell us the net effects? When we say inverse of the community matrix, what exactly is a community matrix. To some community ecologists, it is the site-species matrix, which is definitely not what the authors meant here.

- The term "disordered interactions" are coined multiple times, but during the Intro I wasn't left hanging until reading it again in the Discussion. Can it be defined more early on?

- In the last Intro paragraph, what does "dynamical regime of species sorting" mean?

- Does the last summation term (N of y - N of x) imply mass effect?

- Figure 1e is interesting. Do you think it is true / would help to highlight that the dispersal limitation on more superiour competitors is also crucial to the persistence of an inferior focal species in a patch?

- Perhaps refer to Equation 2 in Section 1.4, because I found it confusing what the actual inference is.

Then, throughout the manuscript I constantly thought about multicollinearity, from a data perspective. I am not able to pinpoint this to any place of the manuscript, but let's say the 2nd Discussion paragraph. Here the authors mentioned "...if the same environmental conditions predictably lead to the same species composition..." Under environmental tracking, if one sample communities in the field they would need to face strong autocorrelation / multicollinarity between neighbour densities and the abiotic covariates. Is this true in the authors' simulated data too? Is this the statistical phenomenon that leads to confounding variables / non-indentifiability during the inference process? These are the questions that haunts my sleep whenever I fit a model including both neighbour densities and environment, so I wonder if the authors have any insights. From the top of my head, I also couldn't think of any literature that points out that aggregated statistics may be robust against multicollinearity even when individual coefficients are not, so I wonder if this is another way of highlighting the novel contribution of this work.

The text on equilibrium and "letting indirect effect play out" in the Discussion is really insightful; I hope it was highlighted more clearly / early on.

Thank you.

​​

In the manuscript “Getting More by Asking for Less: Linking Species Interactions to Species Co-Distributions in Metacommunities”, the authors ask if one can infer species interactions from spatial occurrence data.  More specifically, they first create such data by running a metacommunity model with a known species interaction matrix and then try to recover that matrix using regression techniques, as one would do when mining empirical data. Because their model considers the simplest of cases (linear per-capita growth, homogeneous dispersal), the test of inferential capacity is assumed to be generous and comradely: applying more complicated models can only lead to poorer results. Overall, they find that predictions of the whole species interaction matrix are poor, while estimates of its first two moments are more accurate. In general, the manuscript is very well written, the setup is smart, and the results are novel and (as far as I can see) correct.
 
I have the following comments:
 
1.1: “their abundances are set to zero, and only allowed to vary if their net growth rate dlogNi/dt becomes positive”. At first, confusing to me (I guess the verb “vary” got me sidetracked).
 
As I understand the simulation protocol, there is no regional equivalence: some species are inherently more fit than others. Intuitively (but I did not think this through properly), I’d say regional equivalence without dispersal (species do differ at the pixel level, but each species persists in the same fraction of pixels) affects inference of species interactions because it affects to what extent the landscape represents the experiment one would need to infer interactions. 
 
I wonder why species in the dispersal limited case were only seeded in 50% of the patches. I’d think that seeding them everywhere (in combination to regional equivalence) can produce patterns that permit more reliable inference. 
 
I have the same intuition as the authors that more complex models (displaying more complex dynamics) will not improve inference. However, I’m not sure how model complexity influences inferential capacity vs. simply the number of species (assuming the objective is to parameterize each and every specific interaction using field data). For any model, the number of parameters to estimate would scale exponentially with the number of species, and so one needs a larger design to identify them. If there is a lower limit to inferential capacity, one could ask if inference in large systems will hit rock bottom just the same in complex as in simple models. I’m not suggesting this issue be addressed in this paper, I’m just proposing the authors to ponder this thought.
 
The direct and indirect effects are interesting and I think identical to the definitions of Zelnik et al (https://doi.org/10.1101/2022.12.29.522189), so it may be useful to refer to that paper in section 1.3.
 
It might be worthwhile to discuss in what practical applications it would it be sufficient to know only the statistical properties of the interaction matrix, and what these first two moments could tell us about the system’s dynamic behavior and longer-term outcome. In other words: does what one gets (by asking for less) address one’s needs?
 
Fig1: very nice figure; maybe try to find a way to highlight that the tiles in the lower panel are results for different species.
 
Fig4: I don’t think the color codes are needed for stdA.