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):
- Intercept only model, where allochthony/TP does not vary with
any predictor.
- ‘Size’ only model, where allochthony/TP is best predicted by
individual size.
- ‘Location’ only model where allochthony /TP differs
between catchments but shows no appreciable change over time:
- Univariate ‘Temporal’ model, where allochthony /TPchanged significantly over time across all catchments.
- ‘Size’ controlled ‘Temporal’ model, where changes inallochthony /TP over time is significant after
controlling for differences in individual body size.
- ‘Location’ controlled ‘Temporal’ model, where changes inallochthony /TP over time is significant after
controlling for catchment identity.
- ‘Size’ and ‘Location’ controlled ‘Temporal’ model, where changes inallochthony /TP over time is significant after
controlling for location and individual body size.
- ‘Size’ interaction model where trends in allochthony /TPchange over time differs between individuals of varying body sizes.
- ‘Location’ interaction model, where trends inallochthony /TP change over time differs across
catchments.
- ‘Location’ controlled ‘Size’ interaction model, where trends inallochthony /TP change over time differs across
individuals of varying body sizes, after controlling for catchment
identity.
- ‘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.