Data were collected from a total of 127 0.1 ha circular temporary sample plots by using a random- systematic network laid out in the field. Specifically, we set up a 0.1 km x 0.1 km grid across the study region and selected 130 grid intersection plots at random. The sample plots were established in sites with no evidence of disturbances to minimise the noise in the response variable, removing three plots from the initial selection. Within each plot, the diameter of all tree species with DBH > 7 cm and the total height of all beech trees were measured by using caliper and Vertex IV (Haglöf, Sweden), respectively.
Since environmental changes to the soil are more strongly reflected in the topsoil, five soil samples of the top 10 cm below the litter were randomly taken within each plot using core soil sampler. The soil samples were then mixed and analyzed in the laboratory. Roots, shoots and stones were separated by hand and discarded and the air-dried soil samples were then sieved at 2 mm mesh size. Soil organic matter was determined using the Walkey-Black method (Allison 1975). Total nitrogen was measured in the laboratory by the Kjeldahl method (Bremner and Mulvaney 1982). The available K was determined by a flame atomic absorption spectrophotometer (AA500F, PG Instruments Ltd, China). The available P was determined by using the Olsen method (Homer and Pratt 1961). The Bouyoucos hydrometer method was used for determining the soil texture (Bouyoucos 1962). Soil pH and bulk density (at air-dried moisture content) were determined using pH-meter and Plaster (1985) method, respectively. Site factors such as altitude, slope percent and aspect were recorded at each sample location. Aspect, as the azimuth measured from true north, was then converted to a topographic radiation index using the following equation TRASP = [1 – cos((π/180)(θ – 30))]/2 (Alavi et al. 2019). Environmental variables collected as a basis for modelling are summarized in Table 1.
Statistical analysis
Collinearity among environmental predictors was tested by hierarchical cluster analysis using squared Spearman correlations with the Hmisc package (Harrell et al., 2018) in the statistical software R (R Core Team 2018). The variables percentage sand, percentage carbon, percentage organic matter and percentage saturation were hence removed from the set of predictors due to their high correlation with altitude, slope, TRASP, clay, silt, nitrogen, phosphorus, K, C-N-ratio, bulk density, pH as they were deemed ecologically least relevant based on the authors' expert knowledge. The linear model with quadratic and first-order interaction terms was simplified using backward stepwise and Bayesian Information Criterion (BIC), which considers both the goodness-of-fit and model complexity and penalizes model complexity more than the AIC does (Burnham and Anderson 2002). The model selection step was carried out on the MLR to identify a minimal adequate model structure for all model types. The residuals of the MLR were normally distributed (for details see supplementary information).
Simultaneous autoregressive (SAR) models assume that the response variable at each location i, conditional on the value of explanatory variables at i, depends on the other response variables at neighboring locations j (Haining 2003). SAR models enhance the linear regression model with an additional term that combines the spatial autocorrelation structure of observations in data. In simultaneous autoregressive models, the neighborhood relationship is formally expressed in an n × n matrix of spatial weights (W), in which elements (wij) represent a measure of the connection between locations i and j. In the present study, we used three different SAR models including the spatial error model (SEM), spatial lagged model (SLM) and spatial Durbin model (SDM = spatial mixed model). The SEM models assume that the autoregressive process is in the error term. This is most probably in cases when spatial autocorrelation is not completely explained by the predictor variables (Diniz-Filho et al. 2003). The SAR spatial error model takes the form
Y = Xβ + λWµ + ε ,
where λ is the spatial autoregression coefficient, W is the spatial weights matrix, β is a vector representing the slopes associated with the explanatory variables in the original predictor matrix X, and ε represents the (spatially) independent errors.
The SLM models suppose that the autoregressive process affects only the response variable. It takes the form
Y = Xβ + ρWY + ε
where ρ is the autoregression parameter, and the remaining terms are as above.
If spatial autocorrelation can affect both response and explanatory variables, SDM takes the form (Anselin and Griffith 1988):
Y = Xβ + ρWY + WXγ + ε
Here a new term (WXγ) appears in the model, which represents the autoregression coefficient (γ) of the spatially lagged explanatory variables (WX).
In all SAR models, the neighborhood needs to be provided as input, which we consider as a room for arbitrary decisions (as discussed in Bauman et al. 2018 for eigenvector approaches), particularly compared to the generalized least squares model (GLS). GLS is, in principle, just the name for the algorithm used for fitting regression models with pre-specified error covariances and hence also used to fit the SAR models (Pinheiro and Bates 2000). Typically it is referred to, however, as an approach where the spatial covariance structure is usually modelled assuming a simple distance-decay, e.g. exponential or linear (Dormann et al., 2007).
All approaches, MLR, SAR and GLS, were fitted with the free software R. The three SAR model types (SEM, SLM, and SDM) are implemented in the 'spdep' package (Bivand 2011). Determine the neighborhood distances is an important issue for using the SAR function. After that spatial weight matrix calculated by weighting the neighbors with a certain coding scheme. For specifying the best neighborhood distances, we looped through 20–50 m distance settings and the best neighborhood distance was selected as the one with the lowest AIC. The final Simultaneous Autoregressive Regressions were run with a spatial weights matrix based on a neighbourhood distance of 150 m and a row standardized coding scheme ‘W’. For the GLS model, three different correlation structures (corExp, corGaus, and corSpher) were specified using the nlme package (Pinheiro et al. 2019).
In order to test the spatial autocorrelation in the model residuals, we employed Moran’s I, which quantifies the variance among residuals as a function of geographical distance. Values of Moran’s I vary between − 1.0 and 1.0, and positive values show observations within a certain distance have a tendency to be similar, negative values indicate dissimilarity, and approximately zero means arranged randomly and independently over space (Bao 2001). To evaluate patterns of spatial autocorrelation in the residuals of all regression models we used Moran’s I assessed by a test statistic (the Moran’s I standard deviate) (Dormann et al., 2007). Spatial correlograms for the mean height of oriental beech trees were prepared based on a neighbourhood from 0 to 150 m. For comparison of the different modeling approaches, we computed different goodness-of-fit measures: (1) percentage of variance explained, as adjusted R2; (2) AIC; (3) root mean square error (RMSE) and of course Moran's I.