In forestry, spatial heterogeneity is theorized as one of the major drivers of biological diversity (Wiens, 1976). Spatial heterogeneity results from the spatial interactions between a number of biotic and abiotic factors and the differential responses of organisms to these factors (Milne, 1991). It may have significant influences on many ecosystem processes at multiple spatial scales (Turner, 1989). The spatial heterogeneity of vegetation patterns (i.e., landscape heterogeneity) is a structural property of landscapes (Li & Reynolds, 1994) that can be defined by the complexity and variability of ecological system properties in space. The spatial heterogeneity of a tree variable (height, diameter, etc.) in a forest stand results from the complex historical and environmental mosaic imposed by competition and systematic environmental heterogeneity (Zhang et al., 2004); which implies a spatio-temporal heterogeneity. Forest researchers have realized that the spatial pattern of tree locations strongly affects (1) competition among neighbouring trees, (2) size variability and distribution, (3) growth and mortality and (4) crown structure (Zhang et al., 2004). Ignoring spatial heterogeneity in forest modelling causes biased parameter estimates, misleading significance tests, and sub-optimal prediction (Anselin & Griffith, 1988).

Global and local regression models have traditionally been used for studying spatially heterogeneous data (Fotheringham et al., 2002). There are important differences between these models, for example, global regression models summarize the characteristics of the spatial pattern over the whole study area and attempt to identify homogeneity (i.e. Ordinary Least Squares). On the other hand, local regression models make explicit the differences in the pattern observed among parts in the study area (i.e. Geographically Weighted Regression) (Fortin & Dale, 2005). This aspect makes global models easy to compute. However, these models are not realistic and the results obtained may not be correct. Local models attempt to identify exceptions and are variable throughout the space. Hence, these models are suitable for specific sample plots but are of limited use for large populations (Fotheringham et al., 2002).

There are several methods for solving the heterogeneity problem, for example, models based on heteroscedasticity specification, generalized additive models, classification and regression trees, linear mixed models (LMM) and geographically weighted regression (GWR). In this document the latter two methods were analysed.

LMM has been widely used for height-diameter modelling (Lappi, 1997; Mehtätalo, 2004; Calama & Montero, 2004). This method estimates fixed and random parameters simultaneously for the same model. The introduction of random parameters into the model, specific to each sampling unit, enables us to model the variability detected for a given phenomenon among different locations, after defining a common fixed functional structure (Lindstrom & Bates, 1990). Mixed models give an unbiased and efficient estimation of the fixed parameters of the model. Furthermore, to correctly apply mixed models in unsampled areas a calibration is required, which is of critical importance when making local predictions (Meng & Huang, 2009).

The second methodology analysed, GWR, is used to explore spatially heterogeneous processes (Brunsdon et al., 1996). The underlying idea of GWR is that parameters may be estimated anywhere in the study area given a dependent variable and a set of one or more independent variables which have been measured at known locations (Charlton & Fotheringham, 2009). Finally, no additional measurement is needed to predict the value of parameters for unsampled locations. GWR extends ordinary least squares regression models by accounting for spatial structure and estimates a separate model for each geographic location within the studied area (Matthews & Yang, 2012). It can be said that GWR is similar to a ‘spatial microscope’ in reference to its ability to measure and visualize variations in relationships, unobservable for non-spatial global models (Matthews & Yang, 2012). GWR was originally used in relation to the economy and civil engineering. In the field of ecology, GWR has been applied to remotely sensed data counts with some interesting results (Wang et al., 2005; Salas et al., 2010; Chen et al., 2012). However, in forest management, the GWR method has predominantly been used for small study areas although it was developed for use in large scale studies.

Tree height is an important variable for forest management in several ways: (1) total height together with diameter at breast height allows us to estimate the total volume of the tree, (2) in regular stands, dominant height together with age are used to estimate the site index and (3) the distribution of tree heights in the stand is very important to describe the vertical structure.

Nevertheless, height measurement can be time consuming and consequently expensive, therefore it is common to use height-diameter models instead. There are two kinds of height-diameter models: local models and generalized models.

Local models are designed specifically for forest stands and are the most precise curves for a particular forest stand (Diéguez-Aranda et al., 2009). However, they have several disadvantages, for example: (1) they require the heights of many trees in the stand to be measured, therefore are very expensive and (2) due to forest heterogeneity (age, quality, density) and the silvicultural state of stands within a given forest, a local function cannot be adjusted to all situations (Diéguez-Aranda et al., 2009). Generalized height-diameter models could improve the results in such forests. These models include stand variables for obtaining precise and unbiased estimations (Gaffrey, 1988). Dominant height is usually included because it is easy to measure and is independent of silvicultural treatments. Other variables such as basal area, density or quadratic mean diameter could also be included (Sloboda et al., 1993; Páscoa, 1987).

The main aim of this paper is to analyse the suitability of GWR for modelling height-diameter relationships in a large area and to compare it with LMM, which is one of the most frequently used methods for including spatial heterogeneity in these models.

To test whether GWR provides as good results as the more frequently used LMM, we applied both methodologies to two different functions: the Petterson (1955) local function and the Cañadas et al., (1999) generalized function, and compared their ability to predict heights when applied to an independent dataset.

Table 3 shows the estimated coefficients for equations 3 and 4, firstly for the LMM and GWR when applied to the Petterson height-diameter function, and secondly when they were applied to the Cañadas I height-diameter function. In this table, the minimum, mean, standard deviation and maximum values of the coefficients obtained for the whole study area (considering the map of coefficients for GWR and the fixed and random effects for LMM) were detailed together with R-squared (R²) and the Akaike information criterion (AIC).

For the Petterson height-diameter function, the coefficient a₀ mean value varied between 0.37 and 0.55 for LMM and also for GWR independently of the bandwidth (between 0.34 and 0.52 for a bandwidth of 1000, between 0.32 and 0.69 for a bandwidth of 500, between 0.34 and 0.71 for a Fixed Kernel bandwidth, and between -0.19 and 0.64 for an adaptive kernel bandwidth). The range of this parameter was lowest when obtained by LMM or by GWR with a bandwidth of 1000 m, but increased for the bandwidth of 500 and the fixed kernel, and reached a maximum for the adaptive kernel. The coefficient a_1, i.e. the slope of the curve, showed higher variation with larger standard deviations and wider ranges than the a₀. So, for LMM the standard deviation of a₁ was 0.48 while for GWR were from 0.65 to 1.39. However, mean values along the study area were similar for all the methodologies, from 1.79 to 2.06. Once again, the largest variations were found for the GWR fitted with fixed or adaptive kernel, proving to be more sensitive to the spatial variability of the data. Moreover, the smallest AIC was obtained with the fixed kernel, representing an optimal bandwidth of 363 m; and with the adaptive kernel, representing an optimal number of trees of 26. However, all GWR methods showed certain incongruent results such as negative or very large parameter values. An example of estimation using the Petterson map of coefficients is presented in Fig. 2, which provides a graphical representation of the previously explained parameter distribution for GWR with a predefined bandwidth of 1000 m.

In the case of the Cañadas-I height-diameter function it can be seen that the coefficient b₀ had slightly larger standard deviations and ranges for GWR than for LMM. Additionally, the mean values were almost the same between both methods. The AIC for this function was smaller when GWR was fitted with a predefined bandwidth of 500 m. For this generalized height-diameter function, the optimal bandwidth calculated by fixed kernel was 755 m, and the optimal number of trees in the adaptive kernel was 65. Again, some incongruent minimum values can be seen for the GWR methodology in Table 3. A representative map of coefficients for the Cañadas I estimation is shown in Fig. 3, and the parameter distribution of GWR is presented with a predefined bandwidth of 1000m.

As response variables in the Petterson and Cañadas-I functions were different, the AIC and R² values were presented in Table 3. However, AIC were not comparable between models.

Tree height is an important variable which is used for estimating stand volume or site quality (Diéguez-Aranda et al., 2009). Forest variables are typically highly heterogeneous and in this paper, two different methods were presented for modelling height-diameter relationships and their spatial heterogeneity. The first method was LMM, in which the height-diameter curves use fixed and random coefficients and which gives very accurate predictions when random effects can be predicted (Nanos et al., 2004). However, to calibrate the LMM prior to use, additional height measurements of trees within the stand are required (Lappi, 1997). The second method was GWR, which is a geostatistical method that can predict height-diameter models without the necessity for additional measurements. This effect is very interesting in case of large inventory areas such as national forest inventories, because the only available data are those collected in the inventory plots. These two methodologies were tested according to a local (Petterson) and a generalized (Cañadas I) height-diameter function as well as with different bandwidths in the case of GWR.

In the case of the Petterson model, it was observed that this method had insufficient predictive quality when GWR or uncalibrated LMM was applied. It is possible that the large heterogeneity of the studied area made it difficult for the Petterson function to reach the necessary quality. This could be due to the fact that the Petterson formula only had diameter at breast height as an independent variable and would appear to be too simple to reproduce the forest heterogeneity. Nevertheless, when Petterson equation was validated with calibrated LMM, a high improvement in the predictive quality was observed. In the case of height-diameter model, random component represents the specific conditions of every stand in the model. It causes better precision of the model for particular location in comparison with the model with fixed parameters. This could be important in a local model such as Petterson, designed specifically for giving precise prediction in individual forest stands (Diéguez-Aranda et al., 2009).

In contrast, the Cañadas I model revealed high stability in the model fitting and high predictive quality in the validation in both LMM and GWR. In these cases, R² values in the validation improve from local to generalized regressions from 0.63 to 0.72 in calibrated LMM. Similarly, MAE and RMSE were lower for the generalized models. The comparison between local and generalized equations is analogous with previous studies which concluded that the inclusion of stand-density variables such as dominant height or dominant diameter in the base height–diameter function increased the accuracy of prediction (Temesgen & Gadow, 2004; Diéguez-Aranda et al., 2009).

The predictive quality of the model was improved when few additional height measurements were used for estimating the random parameters for the validation dataset. Calibration using small subsamples of trees is one of the main utilities of mixed models in forest applications (Calama & Montero, 2004). However, when additional measurements are not available GWR provides predictions with higher quality than the LMM.

In some of the GWR models the presence of odd, extreme values was observed. With the Petterson equation these values appear in every GWR model and did not correspond to specific plots but rather, were scattered randomly in the study area. This could be another indicator of the weakness of these models. However, in the case of Cañadas I, these odd values only appeared in the models with a bandwidth of 500 m and adaptive kernel. Specifically, this odd value was located in plot number 1473, which contained only 3 trees and was located in a relatively isolated area. Odd values have also been observed in similar height-diameter modelling studies (Zhang et al., 2004; Zhang et al., 2008; Zhang et al., 2009). It is important to be aware of the presence of these odd values in order to select the correct model. The inclusion of these values in the interpolation could lead to erroneous results in the GWR coefficient map.

In our study, the observed GWR models gave the same residuals regardless of the calibration used. However, according to Zhang et al., (2009), the GWR with adaptive kernel was expected to produce more desirable model residuals in terms of spatial autocorrelation and heterogeneity than the GWR with fixed kernel model for the plots. There are a wide range of calibration options so it is possible that future modifications of settings could improve the current results. For instance, it would be possible to test different bandwidths, or to add robustness to GWR (Fotheringham et al., 2002), or even to calculate the optimal bandwidth differently or the optimal number of trees with different kinds of weighted functions (Gollini et al., 2015).

The maps of coefficient in Figures 2 and 3 illustrate the spatial heterogeneity of the height-diameter relationship. In both figures it can be seen that there are greater slopes in both models in the upper-left area and in the bottom-right area. Additionally, more flat areas exist in the centre of the studied areas. In Fig. 2, two areas can be observed where a₀ takes the highest values and a₁ the lowest; therefore in Fig. 3, only the left side is visible. It would be interesting to compare these contours with local terrain layers (slopes, kind of soil, vegetation), in order to determine whether there is a relationship between fittings and local conditions. Since the spatial variation of local conditions and competition were taken into account in the GWR model, it showed better fitting results for the height-diameter regression relationship (Zhang et al., 2004).

When comparing the LMM and GWR models in this study, it was found that GWR estimation produces similar errors to the LMM when generalized equation was applied. Analogous results have been found in other studies, supporting our finding that GWR gives similar quality estimations to the LMM methodology and predicts the response variable accurately (Zhang et al., 2005; Zhang et al., 2008, Zhang et al., 2009). However, GWR showed not so good predictive quality with the Petterson function, showing no important improvement over LMM without calibration. We suppose that GWR quality could be improved by more sophisticated calibration, which needs a deeper investigation.

LMM and GWR provide different perspectives and insights for data analysis and model prediction. In LMM, the fixed component explains the impact of different variables as with the ordinary least squares regression (Yang & Huang, 2011), and the random component explains the heterogeneity and randomness by both known and unknown factors (Vonesh & Chiinchilli, 1997). LMM is able to characterize the spatial covariance structures in the data with different geostatistic models (Zhang et al., 2009), and produces accurate predictions for the response variable by accounting for the effects of spatial autocorrelation through the empirical best linear unbiased predictors (Littell et al., 1996; Schabenberger & Pierce, 2002). However, GWR has several very important advantages over LMM, which include the following: (1) It is a useful tool for exploring spatial nonstationarity at different scales (Jetz et al., 2005), which was achieved by changing the bandwidth of study. (2) Its capacity for testing the heterogeneity of the study area (Danlin & Yehua, 2013). The heterogeneity is shown in the coefficient maps (Figures 2 and 3) and can be tested numerically by applying the relative spatial heterogeneity index (Zhang et al., 2009). (3) GWR is able to work with spatial data and every kind of attribute data (Danlin & Yehua, 2013); this is particularly useful in the case of Forest Inventories where attributes of measured trees and their location are known. (4) The residuals for GWR are at least similar as for LMM (Zhang et al., 2004; Zhang & Shi, 2004; Zhang et al., 2009), as confirmed in our research. (5) GWR results, unlike LMM model results, are mappable; these maps of coefficients facilitating the interpretation based on spatial context and known characteristics of the study area (Goodchild & Janelle, 2004). The maps of coefficients in our study are shown in Figures 2 and 3. (6) In contrast to the LMM method, the GWR regression does not require any additional measurements or calibrations.

Two different methodologies were used to model local and generalized height-diameter relationships: the geostatistical method GWR and LMM methodology. Additionally, two different height-diameter functions were selected to fit these models: a local function and a generalized function. Finally, the quality of the estimated models was compared through statistical analysis.

As expected, it was observed that in general, both LMM and GWR provide better prediction quality when applied to a generalized height-diameter equation. However, error analysis revealed that there were no great differences between these two methodologies, evidencing that GWR applied to a generalized function provides results which are at least as good as those obtained using the calibrated LMM methodology. However, GWR has a number of advantages over LMM. For instance, GWR provides information on a continuous map of regression coefficients within the studied area. Consequently, is possible to calculate tree heights in every single location without the need for additional measurements. Furthermore, the GWR method can help us to understand spatial processes and may have interesting applications in the area of forest yield and production.