2) Analyses based on random sampling
In the above-mentioned analyses, the sample size differed markedly among
grain sizes [the number of 400, 800, 1200 m2 and
2500 m2 (sub)plots was 61, 30, 20 and 10,
respectively]. In case that the difference in sample size would affect
the statistical results, we randomly extracted a same number (14) of
subplots for the 400, 800 and 1200 m2 grain sizes,
which were fitted with mixed-effect model and then submitted to model
selection (as mentioned above). This procedure was repeated for 200
iterations. We sampled 14 subplots because this was the least sample
size (which was closer to the 10 sample size of 2500
m2 plots) allowed by the MuMin package when we wanted
to use 10 candidate predictors in a mixed effect model. For the 2500
m2 plots, we did not conduct random sampling analyses
because of the limited sample size (10) compared with number of
predictors (see above), and the results were listed only for reference.
In each iteration, for the variables retained in the model, we used
variation partitioning analyses to decompose the total variations in a
response variable into four components (Borcard et al. 1992): the
pure effect of stand factors (a ) and diversity (b ),
respectively; the joint effect of stand factors and diversity
(c ), and unexplained variations (u ) (see Fig. S3 for an
illustration). Note that the pure effect of diversity is the explanatory
power when the effect of stand factors has been already accounted for
(Schmid et al. 2009, Wu et al. 2015), and thus it may be argued that the
biodiversity effect may be underestimated by this metric (the situation
is also true for stand factors). Consequently, we also reported the
total effect for both diversity (i.e. a + c ) and stand
factors (b + c ). In some iterations, only stand factors
(or diversity indices) were retained in the model, in this situation
both the total and pure effects were equal to the model
R2. In a next step, we calculated the mean pure and
total effects for diversity (or stand factors) over the 200 iterations,
to measure the relative effect of diversity (or stand factors) on each
response variable (i.e. biomass, productivity and their temporal
stability). We used these results for a more robust evaluation of the
relative importance of diversity vs. stand factors (Question 1),
and how their relative effects changed with plot size (Question 2).
Finally, to examine the relative roles of different diversity dimensions
and components (Question 3), we also calculated variable importance for
the variables retained in the model, for each iteration of the 200
random samplings. We conducted hierarchy partitioning analyses
using the “hier.part” function
in hier.part package in R (Walsh & Nally. 2020). In each iteration, we
recorded the importance value of each diversity or stand factor variable
(when a variable was excluded from the model, it’s importance value was
0), and model R2. In the “hier.part” function,
variable importance is calculated as the percentage of independent
contribution to model R2 by each variable in the model
(Walsh et al. 2004). However, the model R2 in
different iterations might differ markedly, and thus the variable
importance was not well comparable among iterations (for instance, FRic
may contribute to 80% in a model with an R2 of 0.3,
but may contribute to 50% of a model R2 of 0.8).
Consequently, we calculated “variable importance” in each iteration as
the product of model R2 and the importance value
obtained by hier.part. Then the importance of a diversity dimension (or
a component) was calculated by averaging the importance of corresponding
diversity indices (e.g. importance of the richness component = averaged
“variable importance” for species richness, FRic and PSR). In a final
step, we averaged the resulting importance values across the 200
iterations to evaluate the relative effects of different diversity
dimensions (or components) on ecosystem functions and stability.
All statistical analyses were performed with R 3.6.0 (R Development Core
Team, 2019).