SuppMat 2: Effects of Network Reconstruction on Food-Web Structure

Methods

Effect of Body Size Sampling Method on Network Metrics

To assess whether the choice of body size sampling distribution (uniform, lognormal, or truncated lognormal) influences the estimated structure of ecological networks, we computed partial eta-squared (η²) values for each network metric within each reconstruction model. This approach isolates the effect of body size distribution while controlling for model-specific variation.

Across all models and network metrics, the effect of body size sampling method was extremely small. Most η² values were effectively zero (<0.01), indicating that the choice of distribution had negligible influence on metrics such as connectance, generality, vulnerability, and various motif counts. Only a few metrics in the ATN model showed slightly higher effects (η² ≈ 0.17 for number of linear chains), but these remained the exception rather than the rule.

These findings justify using any of the tested body size sampling approaches for our simulations, as the structural conclusions drawn from the networks are robust to this methodological choice. Consequently, analyses of network metrics and comparisons across reconstruction models are not confounded by the specific form of the synthetic body size distribution.

Figure S1. Effect of Body Size Sampling Method on Network Metrics

Metrics where η² was notably higher are labelled directly on the bars. The figure highlights that, for the majority of network properties and models, η² values are extremely low (<0.01), confirming that the body size distribution choice has negligible impact on network structure. A few exceptions appear for the ATN model, but these are limited to specific metrics (e.g., number of linear chains).

Downsampling PFIM networks

To assess the optimal downsampling parameter to use so as to allow some variation in network structure we downsampled PFIM metawebs across the following range: 5.0, 10.0, 20.0, 30.0, 50.0, 100.0. We then calculated the Jaccard similarity of the resulting adjacency matrices. Across the parameter range the similarity of networks tended to remain above 0.6 and a downsampling parameter of 100 had a mean similarity of 0.86. We thus used a downsampling parameter of 100 for all of the PFIM networks for all statistical analyses.

Multivariate analysis of network structure

We quantified food-web structure using a suite of macro-, meso-, and micro-scale network metrics capturing global topology, motif composition, and species-level interaction patterns (Table S1). Differences among reconstruction approaches were assessed using a multivariate analysis of variance (MANOVA), with model identity as a fixed factor and the full set of network metrics as response variables. Pillai’s trace was used to assess overall multivariate significance due to its robustness to violations of multivariate normality. To identify the multivariate axes driving differences among models, we performed canonical discriminant analysis (CDA) on the MANOVA model. Canonical variates represent orthogonal linear combinations of network metrics that maximise separation among reconstruction approaches. The contribution of individual metrics to each canonical variate was quantified using canonical structure coefficients (correlations between original metrics and canonical scores).

For visualisation, canonical scores were plotted using linear discriminant analysis (LDA), which yields an equivalent discriminant subspace under equal group priors. Model separation in canonical space was visualised using convex hulls encompassing all network replicates for each reconstruction approach. Univariate analyses of variance and effect sizes (partial η²) were calculated for individual metrics and are reported in the Supplementary Materials for descriptive comparison. Pairwise interaction turnover was quantified using link-based beta diversity, which measures dissimilarity in the identity of trophic interactions between networks and captures differences arising from species turnover or changes in interactions among shared species.

Quantification of extinction simulation outcomes and model concordance

To evaluate how reconstruction framework influenced inferred extinction dynamics, simulated community states were compared against observed or expected reference states using two complementary approaches. First, deviation in continuous network metrics (e.g., connectance, mean trophic level, modularity) was quantified using mean absolute difference (MAD). For each metric and time step, MAD was calculated as:

\(MAD = \frac{1}{n} \sum_{i=1}^{n} | M_{i}^{sim} - M_{i}^{ref}|\)

where \(M_{i}^{sim}\) is the simulated value and \(M_{i}^{ref}\) is the corresponding observed or expected value. MAD was chosen because it provides a scale-preserving measure of deviation that is less sensitive to extreme values than squared-error metrics and allows direct comparison across reconstruction frameworks.

Second, agreement in predicted persistence outcomes was evaluated using a modified True Skill Statistic (TSS) at both node and link levels. At the node level, species were classified as present (persisting) or absent (extinct) in each simulated network and compared to their presence–absence status in the reference community.

At the link level, each possible species pair was classified according to the presence or absence of a trophic interaction in the simulated versus reference network. Thus, link-level evaluation quantified agreement in the retention or loss of specific trophic interactions, independent of overall species richness.

For both node- and link-level classifications, outcomes were assigned as true positives (TP), true negatives (TN), false positives (FP), or false negatives (FN), and TSS was calculated as:

TSS = Sensitivity + Specificity − 1

where Sensitivity = TP/(TP + FN) and Specificity = TN/(TN + FP). TSS ranges from −1 to 1, with 1 indicating perfect agreement, 0 indicating performance no better than random expectation, and negative values indicating systematic disagreement. Because TSS is prevalence-independent, it is appropriate for extinction simulations in which class imbalance may occur (e.g., many persisting species or many absent links).

Results

Effects of Network Reconstruction on Food-Web Structure

Table S1. Descriptive statistics (mean ± standard deviation) of network metrics by model

In [1]:

library(knitr)
library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.2.0     ✔ readr     2.1.6
✔ forcats   1.0.1     ✔ stringr   1.6.0
✔ ggplot2   4.0.2     ✔ tibble    3.3.1
✔ lubridate 1.9.5     ✔ tidyr     1.3.2
✔ purrr     1.2.1     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

readr::read_csv("tables/Table_S1_descriptive_stats.csv") %>%
kable()

Rows: 13 Columns: 7
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (7): statistic, ADBM, ATN, Body-size ratio, Niche, PFIM, Random

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

statistic	ADBM	ATN	Body-size ratio	Niche	PFIM	Random
Complexity	0.779 ±1	0.794 ±1	0.865 ±1	0.895 ±1	0.871 ±1	0.744 ±1
Connectance	0.298 ±0	0.225 ±0	0.172 ±0	0.118 ±0	0.119 ±0	0.221 ±0
Diameter	1 ±1	1.38 ±1	4.37 ±4.4	3.8 ±3.8	3.3 ±3.3	2.66 ±2.7
Distance	2 ±2	2.35 ±2	2.81 ±2.81	2.19 ±2.2	2.5 ±2.5	1.02 ±1
Generality	0.902 ±1	1.03 ±1	0.637 ±1	1.12 ±1	1.72 ±2	0.326 ±0
Max trophic level	3.58 ±4	3.2 ±3	6.25 ±6.2	4.03 ±4	2.48 ±2	5.89 ±5.89
No. of apparent competition motifs	1.39 ±1	1.17 ±1	0.501 ±1	0.323 ±0	0.523 ±1	0.0297 ±0
No. of direct competition motifs	1.45 ±1.5	1.34 ±1.3	0.529 ±1	0.143 ±0	0.0644 ±0	1.62 ±1.6
No. of linear chains	0 ±0	0.0291 ±0	0.749 ±1	0.192 ±0	0.104 ±0	0.0108 ±0
No. of omnivory motifs	1.19 ±1	0.456 ±0	0.141 ±0	0.109 ±0	0.139 ±0	0.181 ±0
Redundancy	9.44 ±9.438	6.86 ±6.857	5.07 ±5.07	3.02 ±3.02	2.92 ±2.9	6.77 ±6.765
Richness	34.8 ±34.75	34.8 ±34.75	34.8 ±34.75	34 ±33.99	34.8 ±34.75	34.7 ±34.71
Vulnerability	0.889 ±1	1.07 ±1	0.69 ±1	0.648 ±1	0.633 ±1	1.7 ±2

Table S2. Canonical discriminant analysis

In [2]:

readr::read_csv("tables/canonical_loadings.csv") %>%
kable()

Rows: 8 Columns: 4
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): Metric
dbl (3): Can1, Can2, Can3

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Metric	Can1	Can2	Can3
connectance	-0.75	0.36	-0.45
trophic_level	-0.31	-0.75	-0.25
generality	0.73	0.58	0.32
vulnerability	-0.86	-0.20	0.38
S1	0.38	-0.60	-0.48
S2	-0.35	0.62	-0.53
S4	-0.81	0.14	-0.11
S5	-0.12	0.74	-0.52

Figure S1. Canonical Loadings

Canonical loadings for the first two canonical variates (CV1, CV2) from the canonical discriminant analysis of network metrics. Arrows indicate the contribution of each metric to the multivariate separation among reconstruction models. Colours denote the scale of each metric: Macro (light brown), Meso (brown), Micro (sienna). Metric labels are shown for the most influential variables.

PERMANOVA Variance Partitioning

To quantify the relative contributions of reconstruction framework and temporal turnover to variation in inferred network structure, we conducted permutational multivariate analysis of variance (PERMANOVA). Euclidean distance matrices were calculated from standardised (z-transformed) network metrics. Reconstruction framework (‘model’) and extinction phase (‘time’) were analysed separately to estimate their total contributions to variance, and in combination to assess interaction effects. Significance was assessed using 999 permutations.

Robustness of model effects after temporal centering

To determine whether the dominance of reconstruction framework reflected absolute structural shifts among extinction phases, we repeated the analysis after centering network metrics within each time bin. This procedure removes mean temporal differences while preserving within-phase structural variation. Even after temporal centering, reconstruction framework explained 84.8% of multivariate variance (R² = 0.848, p < 0.001), exceeding the variance explained in the uncentered analysis. Thus, the strong influence of model identity is not attributable to temporal mean differences, but reflects intrinsic structural divergence among reconstruction frameworks.

Statistical Drivers of Network Variation

Statistical Robustness and Assumptions

Factorial ANOVA assumptions were validated via residual analysis. Despite significant heteroscedasticity (Levene’s test, p<0.001), the perfectly balanced design (n = 100 per cell) and large sample size (N = 2400) ensure the robustness of the F-test. Visual inspection of Q-Q plots and Residuals-vs-Fitted plots confirmed that the distributions were sufficiently symmetric for parametric analysis.

Figure S2. Temporal Trajectories of Network Structure by Model

Detailed shifts in network properties across the four extinction phases, categorised by organisational scale (Macro, Meso, Micro). Each line represents the mean value for a specific reconstruction framework, with error bars denoting standard error. This figure illustrates the “baseline” differences between models—such as the Niche model’s tendency to over-estimate motif counts—and their divergent responses to species loss.

Table S3. Variance Partitioning of Framework, Time, and Interaction Effects.

Summary of the two-way factorial ANOVA results for all eight metrics. Values represent partial eta-squared (\(\eta^{2}_{p}\)), which quantifies the proportion of variance explained by each factor. The dominance of the ‘Model’ term across all scales confirms that framework choice is the primary determinant of network topology.

In [3]:

readr::read_csv("../tables/ANOVA_Results.csv") %>%
kable()

Rows: 8 Columns: 4
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): Metric
dbl (3): Model, Time, Interaction

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Metric	Model	Time	Interaction
Connectance	0.871	0.029	0.171
Max trophic level	0.821	0.322	0.324
Generality	0.959	0.230	0.405
Vulnerability	0.900	0.010	0.168
No. of linear chains	0.959	0.514	0.783
No. of omnivory motifs	0.955	0.673	0.782
No. of direct competition motifs	0.973	0.915	0.863
No. of apparent competition motifs	0.959	0.771	0.605

Figure S3. Model Disagreement (CV%) Across Extinction Phases

Trends in inter-model disagreement, quantified as the Coefficient of Variation (CV%) between framework means. The Y-axis is standardised across panels to facilitate comparison between organisational scales. A characteristic “dip” at the ‘during’ phase in several meso-scale metrics illustrates the structural canalisation effect, where severe species loss forces a temporary convergence in model predictions.

Table S4. Percentage Disagreement Between Frameworks Across Extinction Phases.

Calculated inter-model CV% for each metric at each time step. These data points underpin the bubble sizes in Figure 2 and the trajectories in Figure S3. Note the reduction in CV% for linear chains and omnivory motifs during the peak extinction phase (During).

In [4]:

readr::read_csv("../tables/Model_Agreement_CV.csv") %>%
kable()

Rows: 8 Columns: 6
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (2): statistic, level
dbl (4): Pre-extinction, During extinction, Early extinction, Late extinction

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

statistic	level	Pre-extinction	During extinction	Early extinction	Late extinction
Connectance	Macro	39.68855	32.15691	34.79191	40.54695
Generality	Micro	58.51606	41.01746	48.53126	51.75989
Max trophic level	Macro	41.65860	31.46020	32.30812	39.32523
No. of apparent competition motifs	Meso	74.26024	83.35463	84.28280	79.62819
No. of direct competition motifs	Meso	80.49528	81.63039	79.94947	82.73417
No. of linear chains	Meso	173.66964	131.61905	156.24100	160.33886
No. of omnivory motifs	Meso	119.22142	102.05423	101.35425	123.04302
Vulnerability	Micro	47.92763	41.32709	39.99020	47.67690

Figure S4. Pairwise Framework Comparisons (Tukey HSD)

Heatmap showing significant differences between specific pairs of reconstruction frameworks across each extinction phase. Colours represent the magnitude and direction of the difference (estimate); asterisks (\(*\)) indicate statistical significance (p<0.05). This identifies which specific models drive the high CV values seen in Figure 2.

References