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).