1)Analyses based on all (sub)plots
Prior to data analyses all variables were Z-score standardized to avoid the difference in dimension among variables. To examine Questions 1~3, we explained biomass, productivity and their stability (S_B, S_P) with two groups of variables as follows: 1) Biodiversity, including eight variables depicting the richness, evenness and divergence of species, functional and phylogenetic diversity, as mentioned above; 2) Stand factors, including stand density and DBHmax.
In addition to bivariate analyses, we used model selection based on modified Akaike information criterion (AICc) to obtain the most parsimonious model, so as to identify the major factors affecting ecosystem functions and stability. Model selection was conducted with the “dredge” function in the MuMIn package of R, which select the optimal model based on both the lowest AICc value and the least number of predictors (Bartoń 2016). To evaluate the relative importance for variables retained in the models, we calculated Chi-square values for mixed-effect models (400, 800, and 1200 m2 subplots), while F-value for models of the 2500 m2 plots. The 400, 800 and 1200 m2subplots were split from the 2500 m2 plots, and thus were statistically not fully independent among subplots within a same plot. Consequently, for these subplots we used mixed-effect models with plot as a random effect, which were then submitted to the MuMIn package for model selection. Mixed-effect model analyses were implement with the R package “nlme”. For the 2500 m2 plots, ordinary multiple regression was used. However, because of the limited sample size (10 plots), a maximum of seven predictors were allowed by the model selection procedure with the MuMIn package. As a result, we selected five out of eight diversity indices (i.e. species richness and evenness, PSR, PSV and Fric), which showed stronger correlations with the four response variables, to be used in the initial models of the 2500 m2 plots (together with stand density and DBHmax).