Statistical analyses
We divided our analyses in two parts, each corresponding to one of our research questions. In the first part, we assessed temporal trendsallochthony and TP (henceforth referred to collectively asW ) by fitting 11 models to our data set of 96 predatory fish specimens (13 species in 6 genera) (n =96, Supplementary Table 1). Each model describes a competing hypothesis where allochthony andTP were tested separately against combinations of the following predictors: 1) time (i.e., categorical predictor with two levels (before or after ) representing temporal identity); 2)size (i.e., continuous predictor representing the scaled standard length of fish individuals); and 3) catchment (i.e., categorical predictor with four levels representing river catchment identity). In addition to these predictors, we also included taxonomic information (genus) as random intercept and/or random slope terms to control for possible phylogenetic differences in W(allochthony /TP ).
Our competing models and the hypotheses they describe are as follows (details of models in Supplementary Table 2):
  1. Intercept only model, where allochthony/TP does not vary with any predictor.
  2. ‘Size’ only model, where allochthony/TP is best predicted by individual size.
  3. ‘Location’ only model where allochthony /TP differs between catchments but shows no appreciable change over time:
  4. Univariate ‘Temporal’ model, where allochthony /TPchanged significantly over time across all catchments.
  5. ‘Size’ controlled ‘Temporal’ model, where changes inallochthony /TP over time is significant after controlling for differences in individual body size.
  6. ‘Location’ controlled ‘Temporal’ model, where changes inallochthony /TP over time is significant after controlling for catchment identity.
  7. ‘Size’ and ‘Location’ controlled ‘Temporal’ model, where changes inallochthony /TP over time is significant after controlling for location and individual body size.
  8. ‘Size’ interaction model where trends in allochthony /TPchange over time differs between individuals of varying body sizes.
  9. ‘Location’ interaction model, where trends inallochthony /TP change over time differs across catchments.
  10. ‘Location’ controlled ‘Size’ interaction model, where trends inallochthony /TP change over time differs across individuals of varying body sizes, after controlling for catchment identity.
  11. ‘Size’ controlled ‘Location’ interaction model, where trends inallochthony /TP over time differs across catchments, after controlling for individual body size.
All 11 competing models were parameterised with a Bayesian approach on the rjags *4.6 statistical package (Plummer, 2016). For each model, we ran 100,000 iterations (burn-in=10,000) on four parallel chains (thinning=1), with vaguely informative priors for both fixed and random effects. We compared the models using the Widely Applicable Information Criterion (WAIC) and the Efficient Approximate Leave-One-Out (LOO) indices (Vehtari et al., 2017) where lower values indicate greater parsimony.
In the second part of our analyses, we assessed the role of forest loss in driving food web changes. As such, we focused on food web measures (i.e., allochthony and/or TP ) reflecting different temporal trends across the catchments (i.e., if models (ix) or (xi) were best supported by the data in the first part of our analyses). This portion of our analyses required the aggregation of our individual-level response variable, W , because our forest loss metrics were measured at the catchment-basin scale. We calculatedΔ \(\overset{\overline{}}{W}\), which represents the mean pairwise difference in allochthony(Δ \(\overset{\overline{}}{\text{allochthony}}\)) and/or TP (Δ \(\overset{\overline{}}{\text{TP}}\)) between fish individuals of the same species (to control for potential phylogenic confounders) from respective time points. We calculatedΔ \(\overset{\overline{}}{W}\)(n =4) from a total of 141 individual-level, species-specific pairwise comparisons of Wbetween the before and after time point (ΔW ). Considering the significant sample size restrictions, we maximised the rigour of our community-level analyses by allocating 25% percent of the individual-level pairwise data points for model-testing (described later) before the subsequent aggregation of the remaining ‘training’ data subset. The process was randomised and catchment-specific, meaning that we allocated an equal percentage of individual-level data for model training/parameterisation (75%) and model testing (25%) for every catchment. We tested Δ \(\overset{\overline{}}{W}\) against the following forest loss metrics described earlier:ΔFtotal ; ΔFnet ;ΔFratio ; ΔFsub.total ;ΔFsub.net ; andΔFsub.ratio .
In addition to forest loss metrics, we also testedΔ \(\overset{\overline{}}{W}\) against other ‘null’ predictors unrelated to forest loss. These were: 1) Atotal , a continuous variable representing the total area of river catchments; 2) Asub.total , a continuous variable representing the size of the catchment area immediately upstream of our sampling points; and 3) Δ \(\overset{\overline{}}{Q}\), a continuous variable representing the mean pairwise change in pH measurements between both time points.
We had nine community-level models in total (i.e., three ‘null’ models and six predictive models using forest-loss metrics) which can be summarised as follows:
Δ \(\overset{\overline{}}{W_{m}}=\beta_{0}+\left(\beta_{1}\right)X_{m}\)…..(eq. 2),
where β1 represents the fixed-slope term describing the relationship between Δ \(\overset{\overline{}}{W}\) and the predictor variables which we collectively annotate here as X(i.e., Δ \(\overset{\overline{}}{Q}\),Atotal , Asub.total ,ΔFtotal , ΔFnet ,ΔFratio , ΔFsub.total ,ΔFsub.net , andΔFsub.ratio ), for the m -th river catchment. All predictors were scaled to facilitate model convergence. We parameterised the models with a Bayesian approach using the same settings as with our individual-level models. While Bayesian techniques are generally more robust to small sample sizes (Kery, 2010), we sought to improve the accuracy of our findings by repeating our parameterising procedure 100 times. We used a different random subset of individual-level pairwise data (i.e., 75% of total pairwise comparisons) for the calculation of Δ \(\overset{\overline{}}{W}\)by repeating the randomised draws at the start of every iteration. We recorded the following indices for model comparison: 1) mean slope coefficients values (β1) describing the relationship between Δ \(\overset{\overline{}}{W}\)and the predictor variables; 2) WAIC and LOO as measures of model parsimony (Vehtari et al., 2017); and 3) Root Mean Square Error (RMSE) as a measure of model error. We calculate (3) using residuals between predictions made using β1 coefficients (parameterised using the training data subset) and observed values in the testing data subset.