 Research Paper
 Open access
 Published:
Determination of optimal tree height models and calibration designs for Araucaria araucana and Nothofagus pumilio in mixed stands affected to different levels by anthropogenic disturbance in SouthCentral Chile
Annals of Forest Science volume 80, Article number: 18 (2023)
Abstract
Key message
Here, we present a workflow for determining the optimal tree height model and calibration design for forests affected to varying degrees by anthropogenic disturbance. For mixed AraucariaNothofagus forests, tree height predictions in newly surveyed stands are most accurate and effective when the height of up to five random trees is measured to recalibrate predefined nonlinear mixedeffects models.
Context
AraucariaNothofagus forests in Chile are affected by anthropogenic disturbances such as intentional forest fires, grazing, and seed harvesting, causing forest structure to become more heterogeneous. This also challenges tree height predictions, which are required for yield estimations, carbon accounting, and forest management, since height measurements of standing trees are often considered too costly, difficult, and imprecise.
Aims
How does the structure of these forests vary by different levels of anthropogenic disturbance? Which models for estimating tree height of Araucaria araucana and Nothofagus pumilio are most reliable and generally usable? And considering their application in stands they have not been fitted to, which calibration design is optimal for these models?
Methods
Twelve stands were surveyed and classified into four different intensities of anthropogenic disturbance. In 25 to 36 plots per stand, horizontal point sampling measurements of stem diameter as well as of height of selected trees were carried out. Different quantitative standlevel properties were calculated to determine forest structure, which was compared among stands by cluster analysis. To identify the optimal heightdiameter (H–D) model, simple models including diameter only as well as generalized models including stand variables were tested, each additionally extended by a nonlinear mixedeffects (NLME) modeling framework accounting for nested and random effects. To further determine tree height in new stands, the optimal model calibration design was identified involving the empirical best unbiased predictor technique.
Results
Forest structure greatly varied among stands affected by different levels of anthropogenic disturbance, which challenged the development of tree height prediction models. Of all the simple H–D models considered, the Gompertz model was the best for A. araucana and the Näslund model for N. pumilio. The models progressively improved by adding stand variables and using NLME techniques. However, our final model comparisons indicate that a calibrated simple NLME model without stand variables should be preferred. It was further found that the optimal calibration design is to use five randomly selected trees.
Conclusion
Although anthropogenic disturbances can have a complex effect on heightdiameter relationships, the same H–D model can be used for stands representing different anthropogenic disturbance levels and recalibrated by costeffective measurements.
1 Introduction
In predicting future ecological dynamics, an understanding of the influence of anthropogenic disturbance on forests is of crucial importance, as these can weigh considerably in driving tree community dynamics — even compared to climate change (Danneyrolles et al. 2019). Araucaria araucana (Mol.) K. Koch, also known as the monkey puzzle tree, is an endemic tree species in the mountains of SouthCentral Chile and is typically mixed with Nothofagus pumilio (Poepp. & Endl.) Krasser, Nothofagus dombeyi (Mirb.) Oerst, or Nothofagus antarctica (G. Forst.) Oerst. (Veblen 1982). It is currently listed as an endangered species in the IUCN Red List of Threatened Species and was declared a natural monument in 1990, with logging completely prohibited (Fuentes‐Ramirez et al. 2020).
In Chile, there are 254,217 ha of A. araucana natural forests. The forest composition changes with topography (i.e., aspect), altitude, precipitation, natural disturbances, etc. The mixed A. araucanaN. pumilio forests are mainly found between 1000 and 1600 m asl in the Andes Cordillera, but also mixed A. araucanaN. antarctica forests could be found at the same altitude. Around 48% of the A. araucana dominated and natural forests in Chile are protected as stateprotected wild areas (SNASPE) (Hernández et al. 2022). These protected forests remain minimally affected by anthropogenic disturbance, at least in the last 30 years since the protection has been in place. The other 52% of the forest are affected in their structure and forest composition because of different levels of anthropogenic disturbances such as humancaused wildfires, seed harvesting and gathering, decreased regeneration by livestock pressure mainly by cattle and goats, illegal firewood harvesting, land use changes as well as the introduction of invasive plant species (i.e., pines), and seed predation by invasive animal species (González and Veblen 2007; ZamoranoElgueta et al. 2012; Molina et al. 2015; Hernández et al. 2022). Thus, a large proportion of AraucariaNothofagus forests are affected to varying degrees by anthropogenic disturbance. These are often associated with negative impacts on successional processes and lower canopy cover values (Echeverría et al. 2007; Kutchartt et al. 2022). One of the clearest signs of degradation is the lack of natural regeneration (Premoli et al. 2013; FuentesRamírez et al. 2019). Such anthropogenic disturbances can therefore challenge forest functioning and also its management. This implies the need for new research on the structure of forests under different anthropogenic disturbance levels, on management practices, and on strategies and tools for forest inventories, in particular for AraucariaNothofagus forests.
In this article, we are examining and comparing common and more advanced options to develop tree height models, which are important tools for forest inventory and planning, using A. araucana and N. pumilio as examples. In Chile, these two species, a slowgrowing relict conifer and a relatively shortlived shadeintolerant hardwood, are both important species with high ecological and economic importance. Technical reports by Gayoso (2013a, 2013b) presented heightdiameter models for these two species, but they have not been assessed by crossvalidation. This issue refers to a general question that also applies to models for other tree species: How well can existing models be applied to stands that have been altered in their structure by anthropogenic disturbances? It is expected that anthropogenic disturbances will cause changes in stand structure and generate variability that will lead to large errors in model predictions if such standspecific differences are ignored.
Structure is the most obvious characteristic of a forest stand and can be obtained by measuring the diameter at breast height (DBH, 1.3 m aboveground), height, or crown radius of trees in the stand. The heightdiameter relationship is mainly used to describe the vertical structure of forest stands (González et al. 2001). DBH can be easily and accurately measured during fieldwork. However, height measurements are often considered difficult and costly because they are timeconsuming, visibility is often impaired, and the probability of human measurement errors is high (Colbert et al. 2002). It is therefore a common approach in forestry to measure the DBH of all trees and the height of a certain number of sample trees and then construct a heightdiameter model (H–D model) and use it to predict the missing height of trees (Adame et al. 2008) to obtain the input values needed for subsequent studies.
Different models have been proposed to determine tree heightdiameter relationships for different tree species and regions. Curtis (1967), for example, summarized many available linear H–D models. Because of the theoretical maturity of linear models at that time, they were often the preferred choice. However, today, nonlinear models are widely used to predict tree height because of the nonlinear nature of the heightdiameter relationship and modern statistical software can easily fit such models providing accurate predictions (Huang et al. 1992; Paulo et al. 2011). Metrics such as the crossvalidationbased rootmeansquare error and the percentage relative standard error have been reported to be very useful for selecting accurate and robust models with different purposes (Sileshi 2014; Kutchartt et al. 2021).
In addition, H–D models can be described as simple (or local) and generalized (or regional) (Gollob et al. 2018). Simple H–D models, which interpret height as a function of DBH only, are easily applicable without additional measurement, but do not consider the variability of heightdiameter relationships among stands. However, heightdiameter relationships vary between stands and are influenced by the growing environment and stand conditions. Even within the same stand, heightdiameter relationships can change over time (Schmidt et al. 2011). Generalized H–D models use stand variables to estimate local variation of the heightdiameter relationship, avoiding to fitting separate simple H–D model for each stand (Ciceu et al. 2020). Commonly used stand variables are quadratic mean diameter, basal area per hectare, and stand density, which require no additional measurements beyond the full DBH of the stand, making the generalized and simple models equally applicable (Mehtätalo et al. 2015). Calama and Montero (2004) summarized four approaches to include stand variables into models, of which the twostage approach (Ferguson and Leech 1978) has often been used in later studies (e.g., Ciceu et al. 2020). The twostage approach involves fitting a separate heightdiameter curve for each sampling unit (plot or stand) in the first stage to obtain estimates of the parameters and then using the stand variables as covariates to explain the parameters in the second stage. An alternative to such a stepwise approach is to do an exhaustive search across all combinations of stand variables.
The ordinary least square (OLS) technique was the first tool for modeling heightdiameter relationship (Huang et al. 2000). However, H–D models developed with this method usually faced some problems. For example, the data used were obtained from longitudinal measurements or were spatially or temporally correlated, violating the assumption of random and independent observations for modeling and thus leading to biased estimates of parameter confidence intervals (e.g., Dorado et al. 2006; Özçelik et al. 2018; Ciceu et al. 2020; Ogana 2021). To deal with such autocorrelation problem from data, the nonlinear mixedeffects (NLME) models usually fitted by maximum likelihood (ML) method were used (Ercanlı 2020). The parameters of NLME models are divided into two groups in its model structure: fixed effects and random effects parameters (Pinheiro and Bates 2000). The fixed effects parameters show trends in height common to the stand in general, while the random effects parameters account for differences among stands and define variation in the heightdiameter relationship (Calama and Montero 2004). And NLME technique can be applied both to simple and generalized H–D models (GómezGarcía et al. 2015). Furthermore, if height and DBH measurements are taken in a new stand, the random parameters of the mixed effects model can be easily calibrated for that given stand (Mehtätalo et al. 2015). All these model types could be helpful in addressing the variability caused by anthropogenic disturbances at the stand scale. For estimating height for A. araucana and N. pumilio in SouthCentral Chile, no such NLME models have been developed so far.
We hypothesize that forest structural characteristics are affected by the level of anthropogenic disturbance and that this has to be considered in the development of tree height models, which is exemplified here for A. araucana and N. pumilio in SouthCentral Chile. In order to facilitate model application in new stands, our objectives are to use independent data to assess the prediction performance of varying H–D models and to determine their best calibration design.
2 Material and methods
2.1 Study area
The data were collected from 12 stands distributed in the mixed AraucariaNothofagus forests of the Andes Cordillera in the Araucaria region (Fig. 1). The study area ranged in elevation from 1304 m (stand 5) to 1691 m (stand 9) above sea level. The slope in the stands was gentle or flat, and the tree layer was mainly composed of A. araucana and N. pumilio. N. dombeyi was also present, but it occurred rarely and in only two of the stands. Based on CR2MET Version 2.0 gridded climate data, the average annual precipitation ranged from 1930 mm (stand 7–12) to 2582 mm (stand 6) between 1979 and 2019. Over the same time interval, the mean maximum summer temperature ranged from 18.1 °C (stand 6) to 21.5 °C (stand 7–12), and the mean minimum winter temperatures ranged from 5.8 °C (stands 2 and 4) to 6.7 °C (stands 7–12) (Boisier et al. 2018). The dominant soil type in the study area is andosols (IUSS Working Group 2015), which are formed in volcanic tephra and are typically quite young and fertile. CIREN (2002) reported more detailed information for the soils in the western part of the study area. They are characterized by a sandyloamy texture, moderately rapid permeability, and excessive drainage. While moderately acidic (pH from 5.9 to 6.1) and rich in organic carbon in the top layer, the effective cationexchange capacity is very low and the base saturation only 2–5%.
Araucaria araucana (Molina) K. Koch is often referred to as the monkey puzzle tree because of its extremely straight, cylindrical bole and whorled branches (Veblen 1982). It is a longlived, slowgrowing but huge relict conifer native to SouthCentral Chile and southwestern Argentina (Mundo et al. 2013). A. araucana can reach a height of 45 m, a diameter of up to 2 m, and a maximum age of at least 1300 years (Montaldo 1974; Premoli et al. 2013). It is characterized by a thick insulating bark that can exceed 15 cm and a concentration of leaves in the crown, often more than 15 m above the ground, which makes this species resistant to fires (Dickson et al. 2021) which have been shown to be a key disturbance controlling the dynamics of A. araucanaNothofagus forests (Veblen 1982).
Nothofagus pumilio (Poepp. & Endl.) Krasser is a broadleaved tree native to Chile and Argentina, commonly known as lenga beech in English (Barstow et al. 2017). This species can reach a height of over 30 m, a diameter of 1.7 m, and an age of 350 years and is a highly valuable source of quality timber for the timber industry (Magnin et al. 2021). Overall, N. pumilio appears to be less healthy and vigorous under shaded conditions compared to other Nothofagus species such as N. betuloides (Rebertus and Veblen 1993). The relatively short life span and shade intolerant characteristic make this species an obligate seeder that will regenerate rapidly from seed postfire (Dickson et al. 2021).
The 12 selected stands fell into four defined levels of anthropogenic disturbance that were in detail described and assessed by Hernández et al. (2022): none (stands 4, 5, and 6), low (stands 7, 8, and 9), medium (stands 10, 11, and 12), and high (stands 1, 2, and 3). Classified as “none” were protected stands such as those in the Conguillio National Park and the Malalcahuello National Reserve, which remained minimally affected by anthropogenic disturbance at least in the last 30 years. The “low” level was given for stands affected by seed extraction, cattle ranching, and old logging operations. Classified as “medium” were stands that also faced seed extraction in combination with browsing, mainly of cattle, and active logging of N. pumilio. The “high” level was assigned to stands disturbed by intensive grazing and seed collection. These were also affected by a large fire in 2002 (17 years prior to data collection) (González and Veblen 2007), which adds a higher effect of natural disturbance that is distorting the anthropogenic influence there.
2.2 Data collection
The selected stands were systematically sampled using a randomly selected point in each stand as a starting plot from which the entire stand was covered with plots 30 m apart. The number of plots per stand ranged from 25 to 36. At each plot, horizontal point sampling was carried out using a basal area factor of 4 m^{2} per hectare. DBH was measured for all selected trees, and height was measured for onethird of them using a Haga device. All heights mentioned in this study are the total tree heights.
The complete database involves three tree species, A. araucana, N. pumilio, and N. dombeyi, and contains valid data for a total of 3124 trees, of which a total of 873 trees were measured for both DBH and height (Table 1). All these data were considered in the calculation of the stand variables. Of these, 451 heightdiameter data pairs of A. araucana and 380 data pairs of N. pumilio were extracted for the study of the heightdiameter relationship for these two species. In order to develop a reliable heightdiameter model for each species, data from stands 2 and 12 were used as a validating database, and data from the remaining 10 stands were used as a fitting database (Fig. 2). For N. dombeyi, not enough data pairs were given. The complete database has been published on Zenodo (Zhou et al. 2022).
2.3 Heightdiameter model development
The definition of appropriate height diameter models can be a challenging task. Here, we have taken up recommendations by Tischer et al. (2020). We begin with an exploration of simple models and continue by examining model extensions that account for covariates and random effects (Fig. 3). First, the bestperforming simple H–D model in the fitting database was selected for each of the two species from 16 alternative simple models (Mehtätalo et al. 2015) that were widely used by common statistical evaluation criteria and residual diagnostics. Then, 25 alternative stand variables were calculated and included as covariates in the simple models by three different approaches. The stand variable or combination of stand variables that brought the most improvement to the fit and prediction was selected, thus extending the simple model upwards to the generalized H–D model. On this basis, the spatial correlation of the four different sites was considered by the nonlinear extra sum of squares (Huang et al. 2000). At last, nonlinear mixedeffects models were developed in three steps according to Fang and Bailey (2001) and Dorado et al. (2006) to address the effects of different stands on the H–D relationship. After development, the models were calibrated using eight different calibration designs to check how well the models performed in predicting tree height in new stands and to find the optimal calibration design.
A description of the specific technical processes for each step can be found in the annex: Detailed workflow for heightdiameter model development. As a result of these processes, a total of four H–D models were developed for each tree species (Table 2). The R script used and an exemplary analysis are provided through a GitHub repository (Zhou and Zwanzig 2022).
3 Results
In the following, first, the results of the kmeans clustering analysis of the different stand variables will be given to see if the four different levels of anthropogenic disturbance are clustered, i.e., have converging forest structural attributes. Next, the results of model development for each of the two tree species will be presented in subsections based on simple model selection, the inclusion of stand variables, the comparison of H–D models for different sites, the development of the NLME model, and finally on its calibration and height prediction.
3.1 Anthropogenic disturbance levels and stand variables
The kmeans clustering analysis with predetermined 4 clusters did not result in the same classification for the 12 stands as given by the 4 anthropogenic disturbance levels, but some associations seem to be strong (Fig. 4). Stands 1 and 2 from “high” level were clustered into one cluster, as were stands 7 and 8 from “low” level. Both clusters are located on opposite ends of the first axis. All stands from “medium” level were clustered into one cluster, together with stand 5 from level “none.” This cluster is located next to the “high”level cluster but differs more strongly to the last cluster according to stand characteristics loaded on the second axis. The last cluster included stands 4 and 6 from level “none,” stand 9 from “low” level, and stand 3 from “high” level. Overall, stands with no to low anthropogenic disturbance showed a stronger variation than stands with a medium to high level. Among the stand variables used for clustering analysis, \(RNNP\), \(RNAA\), \(NAA\), and \(HdNP\) were the most significant for clustering (p < 0.001), indicating that the number of A. araucana and N. pumilio in mixed stands and their proportions, as well as the DBH diversity of N. pumilio, was important for clustering. \({H}_{dom}\), \({H}_{dom}NP\), and \({H}_{dom}AA\), representing the dominant height in the stands, also had a significant effect on clustering (p = 0.003, 0.006, and 0.013), but of the variables representing the dominant diameter in the stand, only \(DB{H}_{dom}NP\) had a significant effect (p = 0.025). The differences of these stand variables that contribute significantly to clustering can be seen in Fig. 9 in Appendix.
3.2 Heightdiameter model development and calibration
3.2.1 Araucaria araucana (Molina) K. Koch
Simple model selection
Of all sixteen simple models tested, SM12 (Gompertz) was chosen as the simple model for A. araucana. Among the alternative models where PRSE did not show excessive uncertainty in the parameter estimates, SM12 showed the lowest RMSE for both fitting and crossvalidation data, the lowest MAE and the highest R^{2} values (Table 3).
The residual diagnostics showed heteroskedasticity (Fig. 10 in Appendix), as expected for the increasing variability of errors. Different values of \(k\) were tried according to the weighting factor \({w}_{i}=1/{DBH}_{i}^{k}\) , and it was found that the most effective improvement in heteroskedasticity was observed from the residual plots when \(k=1\) (Fig. 10c), similar to other studies (e.g., Huang et al. 1992).
Inclusion of stand variables
The comparison of the three different approaches to construct generalized models revealed different pros and cons. The results of the first approach were not satisfactory. Even when five stand variables were added, the fitting statistics of the model was just similar to that of the model developed by adding only one stand variable in the other two approaches. The second approach brought more improvement to the goodness of fit of the models but was subject to insignificant parameters. The third “exhaustive” approach provided the most improvement in the goodness of fit of the model, regardless of the total number of parameters, but at the cost of requiring much more time to obtain the results than the first two approaches. Therefore, ignoring the time cost, we found the exhaustive approach to be the best choice for including stand variables in the H–D model for A. araucana.
This “exhaustive” approach was processed as follows: after specifying the total number of parameters for the model, the stand variable or combination of stand variables that improved the model fits the most could be derived. They were included in the simple model by parameters, and then, these alternative generalized H–D models were refitted in the fitting database, and their goodness of fit and prediction statistic were calculated and compared with SM12 AA (Table 4). The models explained more and more of the variability as more stand variables were included.
The H–D model developed for A. araucana in this study, after considering model complexity, model performance, and biological interpretation of stand variables, was as follows:
where \({H}_{ij}\) and \(DB{H}_{ij}\) are the height and diameter at breast height of \(j\) th tree in \(i\) th stand, \({H}_{dom}A{A}_{i}\) is the average height of the three thickest A. araucana trees in the \(i\) th stand, and \(RangeDBHA{A}_{i}\) is the difference between the diameters at breast height of A. araucana in the \(i\) th stand.
The effects of these two stand variables on the heightdiameter relationship were simulated (Fig. 5). The effect of \({H}_{dom}AA\) on the heightdiameter relationship was significantly greater than that of \(RangeDBHAA\). As \({H}_{dom}AA\) increased, the tree height also increased, and this effect had a greater impact on thicker trees than on thinner ones. Conversely, as \(RangeDBHAA\) increased, tree height decreased, and this effect was more pronounced for thinner trees and probably negligible for thicker trees.
Comparison of H–D models among different sites
Equation (1) was used as a reduced model when comparing heightdiameter relationships among four different sites. The full model extended by the indicator variable method could be written as follows:
Both the Ftest and the LakkisJones test used in this study showed nonsignificant results (Table 5), which suggested that the same reduced model could be used for the four sites without developing separate generalized heightdiameter models for each site. This also indicated that it was reasonable to assume that there was no spatial correlation at site level. The reduced model (Eq. 1) was the final generalized heightdiameter model developed for A. araucana. There was still heteroskedasticity in this model, but this was greatly improved by using the same weighting factor \({w}_{i}=1/DB{H}_{i}\) (Fig. 11 in Appendix).
Nonlinear mixed effects models
First, a nonlinear mixedeffects model was developed on the basis of the simple H–D model SM12 AA. After comparison of alternative definitions, the parameter was considered to be mixed effect, i.e., having both fixed and random effects. The residual diagnostics showed the problem of heteroskedasticity again (Fig. 12 (left) in Appendix). Of the three alternative variance functions, the power variance function with DBH as the base proved to be the most effective in terms of improving heteroskedasticity.
Similarly, the generalized H–D model (Eq. 1) was refitted using the NLME technique. The parameter was considered to have a mixed effect after comparison. Again, the power variance function was used to improve the heteroskedasticity problem of this model (Fig. 13 in Appendix).
By now, all four models developed for A. araucana in the fitting database were completed, and the parameter estimates and goodness of fit of the models can be found in Table 6. The improvement of the generalized H–D model over the simple H–D model was significant, but the improvement of the generalized NLME H–D model over the simple NLME model was negligible. Comparing the two NLME H–D models revealed that for the generalized one, the value of the variance components associated with the random effects (\({\sigma }_{u}^{2}\)) dropped significantly, but the residual withinstand variance (\({\sigma }^{2}\)) did not drop significantly, suggesting the existence of a pattern of variability that cannot be explained by differences among stands.
Calibration of the NLME models and height prediction
The size of calibration designs had impacts on the prediction performance with both NLME H–D models (Table 7, Table 13 in Appendix). The ONLS fit of the generalized H–D model did not converge in both stands of the validation database.
The values of both MAPE and MSPE indicated an increasing prediction performance of the simple NLME H–D model as the number of sample trees increased. The value of MPE illustrated that the calibrated model produced an overall underprediction of height at less than 4 sample trees and an overall overprediction of height starting with a calibration design of 5 sample trees. The prediction bias (%) of the simple NLME H–D model under different calibration designs were calculated separately for the two validating stands (Table 7). There was a clear difference in the prediction bias (%) in the two stands, with significantly better prediction performance in stand 12 than in stand 2. Moreover, the prediction performance of the model increased in each stand as the number of trees sampled increased. Therefore, considering the cost of sampling and, to facilitate comparisons, the calibration results of the simple NLME H–D model of N. pumilio (see results of Sect. Nothofagus pumilio (Poepp. & Endl.) Krasser below), five trees were selected as the calibration design for the simple NLME H–D model.
The value of the MPE derived from the calibration of the generalized NLME H–D model showed fluctuations (Table 13 in Appendix). The other two statistics both got smaller as the number of sample trees increased. The prediction bias (%) in each stand showed similar results to those of the simple NLME model (Table 7). For stand 12, the prediction bias (%) of the simple NLME model could be even smaller than the prediction bias (%) of the generalized NLME model. Considering the cost of sampling and the calibration results of the generalized NLME H–D model of N. pumilio, one tree was selected as the calibration design for the generalized NLME H–D model.
An example was performed for each of the two models in each of the two validating stands (Fig. 6). As can be seen from this example, the curves for the fixedeffects response pattern and the calibrated response pattern were close, but both were different from the curves obtained from the ONLS.
3.2.2 Nothofagus pumilio (Poepp. & Endl.) Krasser
Simple model selection
The twoparameter SM1 (Näslund) with the lowest MAE, RMSE, and RMSE (CV) and the highest R^{2} among the alternative models where PRSE did not show excessive uncertainty in the parameter estimates was finally chosen as the simple model for N. pumilio (Table 8). The residual diagnostics showed heteroskedasticity (Fig. 14a in Appendix), and the most effective improvement in heteroskedasticity was observed from the residual plots when k = 1 (Fig. 14c in Appendix). Therefore k = 1 was followed in this study.
Studentized residuals can be used to identify outliers, which are defined as points that are far from other observations in one, two or ndimensional space (Zwanzig et al. 2020). Three sample trees with studentized residual values larger than 4 were considered as outliers and were finally excluded by the subsequent model development process, following the guidelines for removing “influential observations” that were found to have too much impact on model development and parameterization (Zwanzig et al. 2020).
Inclusion of stand variables
The comparison of the three different approaches to construct generalized models revealed that the third “exhaustive” approach is also preferred for N. pumilio. The choice of stand variables by the first approach caused problems of collinearity in the refitted model. The second approach selected stand variables that improved the goodness of fit of the models and did not lead to collinearity problems, but the third approach produced a similar goodness of fit with the inclusion of fewer stand variables.
This “exhaustive” approach was processed as described before for A. araucana. The stand variable or combination of stand variables that improved the model fit the most was included in the simple model by parameters, and then, these alternative generalized H–D models were refitted in the fitting database (Table 9). \(NlogNP\) was the logarithmic form of the stems number of N. pumilio. The reason for using logarithms is that \(NNP\) did not conform to a normal distribution and was converted to a logarithmic form that conforms to a normal distribution for possible subsequent studies.
The H–D model developed for N. pumilio in this study, after taking into account model complexity, model performance, and biological interpretation of the stand variables, is as follows:
where \({H}_{ij}\) and \(DB{H}_{ij}\) are the height and diameter at breast height of \(j\) th tree in \(i\) th stand, \({{DBH}_{max}}_{i}\) is the maximum DBH in the \(i\) th stand, and \(D{q}_{i}\) is quadratic mean DBH in the \(i\) th stand, both not related to tree species, and \(NlogN{P}_{i}\) is the logtransformed forms of trees per hectare of N. pumilio in the \(i\) th stand.
The effects of these three stand variables on the heightdiameter relationship were simulated (Fig. 7). The effect of \(Dq\) on the heightdiameter relationship was the most pronounced, followed by \(NNP\) and \(DB{H}_{max}\). As \(Dq\) increased, the height decreased, and this effect was more pronounced for thicker trees than for thinner trees. Similarly, height decreased as \(NNP\) increased, and this effect was substantial for heights in stands with 10–310 N. pumilio trees per hectare while slowly decreasing as the number of N. pumilio trees in the stand got higher again. Conversely, height increased with higher \(DB{H}_{max}\), and this effect was greatest for tree heights with DBH around 50 cm in the stand.
Comparison of H–D models among different sites
Equation 3 was used as the reduced model, and the full model was extended by the indicator variable method. Both the Ftest and the LakkisJones test showed nonsignificant results (Table 5). The reduced model (Eq. 3) was therefore accepted as the generalized H–D model for N. pumilio. The weighting factor \({w}_{i}=1/DB{H}_{i}\) also needed to be added to improve the heteroskedasticity problem of the model (Fig. 15 in Appendix).
Nonlinear mixedeffects models
The parameter \(a\) of the simple H–D model and the parameter \(a0\) of the generalized H–D model were considered to be mixed effect after comparison. The same power variance function was used to improve the heteroskedasticity problem of both the simple and generalized NLME H–D models (Figs. 16 and 17 in Appendix).
By now, all four models developed for N. pumilio in the fitting database were completed. The parameter estimates and goodness of fit of the models can be found in Table 10. The RMSE decreased from the simple to the more complex models, but the overall differences were marginal. The value of the components associated with the random effects (\({\sigma }_{u}^{2}\)) was very similar to zero in the simple NLME H–D model and could be seen as zero in the generalized NLME H–D model. The residual withinstand variance (\({\sigma }^{2}\)) was even getting higher. This suggested that the variability in the heightdiameter relationship in the fitting database of N. pumilio was influenced almost not by differences among stands but rather by differences within stands, i.e., between plots.
Calibration of the NLME models and height prediction
Different calibration designs also had impacts on the prediction performance with both NLME H–D models of N. pumilio (Table 11, Table 14 in Appendix). And the generalized H–D model also failed to converge in the fitting to both stands of the validating database. It was worth noting that the generalized NLME model showed worse prediction accuracy than the simple NLME model. This was also illustrated by their prediction bias (%) in stands 2 and 12 (Table 11). It was reasonable to select five sample trees for the simple NLME model and one sample tree for the generalized NLME model for calibration.
To determine the prediction performance of the NLME models calibrated by the selected calibration design, an example was also performed for each of the two models in each of the two validating stands (Fig. 8). Unfortunately, in these randomly selected examples, the a priori height data did not result in a significant improvement in the height predictions for either model.
Results of model development
Except for the generalized NLME H–D model, which had a large bias (bias > 20%) in prediction of the tree height of N. pumilio in stand 12, the prediction bias largely fell within \(\pm 20\text{\%}\), indicating that the calibrated models performed well (Sharma et al. 2019) in the validating database. Thus, all models were fitted again in the complete database (fitting plus validating), and the resulting model parameter estimates and variance components (see Table 15 and Table 16 in Appendix) can be used to predict tree height for A. araucana and N. pumilio in SouthCentral Chile. The most recommended of these are both simple NLME H–D models.
4 Discussion
This study investigated the effect of different anthropogenic disturbance levels on forest structure and developed suitable H–D models for predicting the height of the two tree species A. araucana and N. pumilio in mixed stands in SouthCentral Chile. As a result of the presented workflow for the development of heightdiameter models, a total of four H–D models were developed for each of the two tree species. These structurally and technically different models vary in their requirements for model input and quality of model output, which are explained and compared in more detail in the following section.
4.1 Forest structure under varying levels of anthropogenic disturbance
The cluster analysis of forest structural characteristics revealed that stands with a medium to high level of anthropogenic disturbance were more similar to each other than stands with no to low impact. This indicates that anthropogenic disturbances have a strong influence on stand characteristics, but when this influence is small or absent, the effects of natural variation, specific site conditions, or history appear to allow greater variation in forest structural characteristics between stands. The variation between plots of the same stand, however, is greater for stands that are more heavily impacted by anthropogenic disturbance, as it has also been shown by canopy scope measurements in these stands, revealing a canopy cover of 87% for stands with almost no anthropogenic disturbance, of 58% for low and 42% for the medium as well as the high level (Kutchartt et al. 2022). The latter stands show signs of forest degradation, known in SouthCentral Chile from unsustainable fuelwood use, grazing, and cutting of the most viable trees of commercially valuable species (Donoso et al. 2022). Various silvicultural techniques exist to rehabilitate such stands, including those aimed at supporting regeneration processes when systems are in a stage of arrested succession (Donoso et al. 2022).
On the other hand, low levels of anthropogenic disturbance may increase forest structural diversity in ways that are reported to be associated with increases in forest productivity, although results for this relationship are quite mixed (Dănescu et al. 2016). For example, thinning is known to be a key silvicultural intervention to improve growth and quality of secondary forests dominated by Nothofagus (Donoso et al. 2022). This may unintentionally apply here to stands with low anthropogenic disturbance.
As forest structure can change significantly at small scales due to undocumented human actions rather than to largescale environmental factors, model predictions for dependent treelevel properties such as the H–D relationship should ideally be based on local optimizations, i.e., calibrated models.
4.2 Model selection
To select the most appropriate function to predict tree height, 16 nonlinear H–D models that have been frequently used in the past were compared. These included seven twoparameter models and nine threeparameter models. Although convergence is considered one of the challenges of fitting 3parameter H–D models (Mehtätalo et al. 2015; Ogana 2021), the nine alternative 3parameter models converged, but their certainty and identifiability for parameter estimation were less than satisfactory (PRSE > 25%). Only the parameter estimates for SM9 (Logistic) and SM12 (Gompertz) were certain and identifiable. The seven twoparameter models did not suffer from this problem. Except for SM5 (Power), their goodness of fit was in fact very similar to that of the threeparameter models.
The simple H–D model chosen for A. araucana was SM12 (Gompertz), a model whose suitability in describing the heightdiameter relationship has been shown in other studies (e.g., Özcelík et al. 2014, Subedi et al. 2018). However, in the study of Zhang (1997), Gompertz was found to underestimate the height of larger trees (BHD > 100 cm). This problem was also observed in this work. All four models for A. araucana underestimated the tree height of large trees with a DBH > 200 cm to varying degrees. The simple H–D model chosen for N. pumilio was SM1 (Näslund), which also provided a satisfactory fit to the majority of the 28 datasets used by Mehtätalo et al. (2015). Our model selection approach demonstrated that combining the PRSE for assessing parameter uncertainty and the crossvalidationbased RMSE for the accuracy in predicting new data facilitates the identification of accurate and robust models. Many of the heightdiameter equations, however, represent very similar functions that are likely to have negligible differences for these and other criteria, as seen here.
4.3 Model improvement
The weighted least squares method that used \({w}_{i}=1/DB{H}_{i}\) as a weighting factor effectively improved the heteroskedasticity problem that occured in the residual diagnostics of both simple models. This can reduce the risk of bias in parameter estimation, but may not substantially improve model performance (Cormier et al. 1992).
One of the most key challenges in developing a generalized H–D model is the selection of additional interpretive stand variables (Raptis et al. 2021). A large range of alternative stand variables involving as well as not involving tree species were prepared here. These alternative stand variables described stand structure, density, species competition in mixed stands, and other important factors that could help modeling heightdiameter relationships.
For A. araucana, \({H}_{dom}AA\) was used as an additional predictor associated with the asymptote coefficient, which was consistent with many previous studies (e.g., GómezGarcía et al. 2015; Raptis et al. 2021). The asymptotic development of the H–D curve is one of the important elements in characterizing the development of the H–D relationship over time, i.e., the asymptotic maximum of the H–D curve increases as a function of age (Eerikäinen 2003). There was a high correlation between age and dominant height, and the inclusion of dominant height in the model also implicitly allowed for age (Dorado et al. 2006). Moreover, although dominant height was reduced to the average height of the three thickest trees in this study, it was still a measure of the maximum height potential of the stand, describing the effect of stand quality on H–D relationship, as reported before in other studies (Calama and Montero 2004; Kershaw et al. 2016). All other factors being equal, tree height increased with increasing \({H}_{dom}AA\), and its effect on the model was greater than that of the other included stand variables. The \(RangeDBHAA\) was included in the simple model as an additional predictor with parameter \(b\), which defines the displacement along the xaxis. According to this, the height of small trees is predicted to be larger for a given DBH, when \(RangeDBHAA\) is lower. This reflects that height growth is typically increased by increased competition in evenaged or evensized stands. When, on the other hand, \(RangeDBHAA\) is larger, small trees may have a lower stress level and tend to exhibit their usual allometric response. The range shows the allometric corridor for the H–D ratio (Pretzsch 2010), i.e., how flexible A. araucana can respond to variable environmental conditions under the influence of different levels of anthropogenic disturbance.
For N. pumilio, a total of three stand variables were included: \(DB{H}_{max}\), \(Dq\), and \(NlogNP\), two of which were related to DBH, representing the maximum and quadratic mean of DBH in the stand, respectively, and did not relate to tree species. In particular, differences in \(Dq\) are affecting height predictions, with larger height predicted in stands with lower \(Dq\). The effect of the density of N. pumilio in the stand was more heterogeneous, with a stronger positive impact on tree height in the range of very low densities. Considering that height growth is rather stimulated by competition and not vice versa, this is counterintuitive, suggesting that the density of N. pumilio reflects and combines the effect of other stand characteristics, all of which are strongly correlated with this variable.
For both tree species, the analysis of the sitebased H–D models showed that it was not reasonable to develop separate models for each site, and that there were no significant differences in the main height growth patterns between the four sites.
The implementation of mixed effects for both the simple and the generalized H–D models improved the goodness of fit of the corresponding models. It is worth noting that although random effects were defined at the stand level in this study, the mixedeffects approach allows for random effects to be defined at different levels depending on the purpose of the study (Bronisz and Mehtätalo 2020), for example, also at the plot level. This was not implemented, as the plots did not have an appropriate sample size.
The NLME H–D models can provide two types of height predictions for new stands not included in the fitting database, a fixedeffects response pattern and a calibrated response pattern. The fixedeffects response pattern of both simple and generalized NLME models presented the lowest prediction accuracy. The generalized NLME model of A. araucana outperformed its simple form, a trend similar to that reported by Raptis et al. (2021), who compared the fixedeffects response pattern of both simple and generalized NLME models in evenaged black pine (Pinus nigra Arn.) natural stands located in Olympus National Park in Greece. However, the two models of N. pumilio gave opposite results. This was because the simple NLME model for N. pumilio showed very low amongstand variance and higher withinstand variance in comparison. This may suggest that the variability in the H–D relationship is not due to differences among stands but between plots, particularly present at high and moderate levels of anthropogenic disturbance, as discussed above. The generalized NLME model not only accentuated this variability among stands by including stand variables but also amplified this variability by the random effects defined at the stand level.
4.4 Model calibration
A calibrated response pattern requires the measurement of heights in the new stands. As the height measurements in our validating database were also a randomly selected sample of stands, we only compared the effect of the different number of trees with random measurements on prediction accuracy. It is clear from our results that as the number of trees measured increased, the mixedeffects model predicted height with increasing accuracy, in line with the results of many studies (e.g., Dorado et al. 2006). We have chosen a calibration design of 5 random trees for the simple NLME models and 1 random tree for the generalized NLME models, which is a combination considering inventory costs and prediction accuracy. However, the number of sample trees to be measured also needs to consider the structure of the new stand. For stands with a homogeneous structure, using a single tree height measurement in the calibration provides high prediction accuracy (Trincado et al. 2007), while for multilayered stands, using a height measurement of at least four trees can result in much lower prediction error (Sharma et al. 2019). In general, the calibration designs for the different tree species given in the different studies mostly have thresholds between 2 and 5 sampled trees (Ogana 2021). This is because the prediction accuracy achieved by increasing the number of sampled trees requires additional inventory costs. It is also worth noting that the accuracy of the calibrated response pattern depends not only on the number of trees sampled but also on their diameter classes (Calama and Montero 2004). Thus, limited by the available database, all that can be considered is the number of randomly selected sample trees, and such a calibration design is likely to be unsatisfactory. In future studies, we can try to explore the best calibration design by selecting some of the thickest, some of the thinnest, and some of the trees close to the mean DBH in the new sample stand and measuring their heights.
5 Conclusion
In conclusion, it is recommended to use the calibrated simple NLME H–D model for A. araucana height prediction. The reason is that considering the inclusion of the stand variable HdomAA in the generalized H–D model, three of the thickest trees need to be measured for the calculation. In contrast, the tree heights of five randomly selected trees result in more accurate predictions for the simple NLME model, and in field work, sample trees can be selected for measurements that are not visually impaired, reducing the difficulty of measuring tree height and omitting increasing inventory costs compared to measuring the heights of the three thickest trees.
For N. pumilio, the use of a calibrated simple NLME model is also recommended, but further research is required to address the sources of the variability in the heightdiameter relationship within each stand, as this will be useful in improving the predictive performance of the model.
The presented methodology of determining the optimal tree height model and calibration design has the potential to inspire future studies aiming to develop tree height models that account for stand variables and mixed effects and that are intended to be calibrated for new stands. This is particularly important given the variable and complex effects of anthropogenic disturbances on stand structure, as observed here for AraucariaNothofagus forests.
Availability of data and materials
The dataset created and analyzed as part of the current study is available on Zenodo: https://doi.org/10.5281/zenodo.7411420.
References
Adame P, del Río M, Cañellas I (2008) A mixed nonlinear height–diameter model for Pyrenean oak (Quercus pyrenaica Willd.). For Ecol Manage 256:88–98. https://doi.org/10.1016/j.foreco.2008.04.006
Barstow M, Baldwin H, Rivers MC (2017) Nothofagus pumilio. The IUCN Red List of Threatened Species 2017.
Bates DM, Watts DG (1988) Nonlinear regression analysis and its applications. Wiley, New York. https://doi.org/10.1002/9780470316757
Boisier JP, AlvarezGarretón C, Cepeda J, Osses A, Vásquez N, Rondanelli R (2018) CR2MET: a highresolution precipitation and temperature dataset for hydroclimatic research in Chile. Geophys Res Abstr 20:EGU201819739
Bronisz K, Mehtätalo L (2020) Mixedeffects generalized height–diameter model for young silver birch stands on postagricultural lands. For Ecol Manage 460:117901. https://doi.org/10.1016/j.foreco.2020.117901
Calama R, Montero G (2004) Interregional nonlinear height–diameter model with random coefficients for stone pine in Spain. Can J for Res 34:150–163. https://doi.org/10.1139/x03199
Ciceu A, GarciaDuro J, Seceleanu I, Badea O (2020) A generalized nonlinear mixedeffects height–diameter model for Norway spruce in mixeduneven aged stands. For Ecol Manage 477:118507. https://doi.org/10.1016/j.foreco.2020.118507
CIREN (2002): Descripciones de suelos materiales y símbolos: Estudio Agrológico IX Región. Publicación CIREN 122, Santiago, Chile, 360 pages. ISBN 956–7153–35–3.
Colbert KC, Larsen DR, Lootens JR (2002) Heightdiameter equations for thirteen midwestern bottomland hardwood species. North J Appl for 19(4):171–176. https://doi.org/10.1093/njaf/19.4.171
Cormier KL, Reich RM, Czaplewski RL, Bechtold WA (1992) Evaluation of weighted regression and sample size in developing a taper model for loblolly pine. For Ecol Manage 53:65–76. https://doi.org/10.1016/03781127(92)900347
Curtis RO (1967) Heightdiameter and heightdiameterage equations for secondgrowth Douglasfir. For Sci 13(4):365–375. https://doi.org/10.1093/forestscience/13.4.365
Dănescu A, Albrecht AT, Bauhus J (2016) Structural diversity promotes productivity of mixed, unevenaged forests in southwestern Germany. Oecologia 182:319–333. https://doi.org/10.1007/s0044201636234
Danneyrolles V, Dupuis S, Fortin G, Leroyer M, de Römer A, Terrail R, Arseneault D (2019) Stronger influence of anthropogenic disturbance than climate change on centuryscale compositional changes in northern forests. Nat Commun 10(1):1–7. https://doi.org/10.1038/s4146701909265z
Dickson B, Fletcher MS, Hall TL, Moreno PI (2021) Centennial and millennialscale dynamics in AraucariaNothofagus forests in the southern Andes. J Biogeogr 48(3):537–547. https://doi.org/10.1111/jbi.14017
Donoso P, Promis A, Loguercio G, Attis Beltrán H, Caselli M, Chauchard L, Cruz G, Peñalba M, Pastur G, Navarro C, Núñez P, SalasEljatib C, Soto D, VásquezGrandón A (2022) Silviculture of South American temperate native forests. NZ J For Sci. 52:2. https://doi.org/10.33494/nzjfs522022x173x
Dorado FC, DieguezAranda U, Anta MB, Rodríguez MS, von Gadow K (2006) A generalized heightdiameter model including random components for radiata pine plantations in northwestern Spain. For Ecol Manage 229(1–3):202–213. https://doi.org/10.1016/j.foreco.2006.04.028
Echeverría C, Newton AC, Lara A, Benayas JMR, Coomes DA (2007) Impacts of forest fragmentation on species composition and forest structure in the temperate landscape of southern Chile. Glob Ecol Biogeogr 16(4):426–439. https://doi.org/10.1111/j.14668238.2007.00311.x
Eerikäinen K (2003) Predicting the height–diameter pattern of planted Pinus kesiya stands in Zambia and Zimbabwe. For Ecol Manage 175(1–3):355–366. https://doi.org/10.1016/S03781127(02)00138X
Ercanlı İ (2020) Innovative deep learning artificial intelligence applications for predicting relationships between individual tree height and diameter at breast height. For Ecosyst 7(1):12. https://doi.org/10.1186/s40663020002263
Fang Z, Bailey RL (2001) Nonlinear mixed effects modeling for slash pine dominant height growth following intensive silvicultural treatments. For Sci 47(3):287–300. https://doi.org/10.1093/forestscience/47.3.287
Ferguson I, Leech J (1978) Generalized least squares estimation of yield functions. For Sci 24:27–42. https://doi.org/10.1093/forestscience/24.1.27
FuentesRamírez A, ArroyoVargas P, Del Fierro A, Pérez F (2019) Postfire response of Araucaria araucana (Molina) K. Koch: Assessment of vegetative resprouting, seed production and germination. Gayana Bot. 76(1):119–122. https://doi.org/10.4067/S071766432019000100119
FuentesRamirez A, SalasEljatib C, González ME, UrrutiaEstrada J, ArroyoVargas P, Santibañez P (2020) Initial response of understorey vegetation and tree regeneration to a mixedseverity fire in oldgrowth AraucariaNothofagus forests. Appl Veg Sci 23(2):210–222. https://doi.org/10.1111/avsc.12479
Gayoso J (2013a). Funciones alométricas para la determinación de existencias de carbono forestal para la especie Araucaria araucana (Molina) K. Koch (ARAUCARIA). Corporación Nacional Forestal. Santiago, Chile. 49 p.
Gayoso J (2013b). Funciones alométricas para la determinación de existencias de carbono forestal para la especie Nothofagus pumilio (Poepp. Et Endl.) Krasser (LENGA). Corporación Nacional Forestal. Santiago, Chile. 39 p.
Gollob C, Ritter T, Vospernik S, Wassermann C, Nothdurft A (2018) A flexible heightdiameter model for tree height imputation on forest inventory sample plots using repeated measures from the past. Forests 9(6):368. https://doi.org/10.3390/f9060368
GómezGarcía E, Fonseca TF, CrecenteCampo F, Almeida LR, DieguezAranda U, Huang S, Marques CP (2015) Heightdiameter models for maritime pine in Portugal: a comparison of basic, generalized and mixedeffects models. iForest 9(1):72–78. https://doi.org/10.3832/ifor1520008
González A, Gabriel J, von Gadow K, Hermosilla PR (2001) Modelización del crecimiento y la evolución de bosques. IUFRO
González ME, Veblen TT (2007) Incendios en bosques de Araucaria araucana y consideraciones ecológicas al madereo de aprovechamiento en áreas recientemente quemadas. Rev Chil Hist Nat 80(2):243–253. https://doi.org/10.4067/S0716078X2007000200009
Hernández J, González V, Promis Á, Corvalán P, Kutchartt E, Pirotti F, Carrer M (2022) Los bosques de AraucariaLenga. Curacautín, Lonquimay y Melipeuco. Alteraciones de hábitat. Universidad de Chile. Andros Ltda., Santiago, Chile. 161 p
Hinkle D E, Wiersma W, Jurs S G (2003) Applied statistics for the behavioral sciences (Vol. 663). Houghton Mifflin College Division.
Huang S, Titus SJ, Wiens DP (1992) Comparison of nonlinear height–diameter functions for major Alberta tree species. Can J for Res 22(9):1297–1304. https://doi.org/10.1139/x92172
Huang S, Price D, Titus SJ (2000) Development of ecoregionbased heightdiameter models for white spruce in boreal forests. For Ecol Manage 129(1–3):125–141. https://doi.org/10.1016/S03781127(99)001516
IUSS Working Group WRB. 2015. World Reference Base for Soil Resources 2014, update 2015. International soil classification system for naming soils and creating legends for soil maps. World Soil Resources Reports No. 106. FAO, Rome.
James G, Witten D, Hastie T, Tibshirani R (2014) An introduction to statistical learning: with applications in R. Springer, New York
Kassambara A, Mundt F (2020) factoextra: extract and visualize the results of multivariate data analyses. R package version 1.0.7. https://CRAN.Rproject.org/package=factoextra
Kershaw Jr, JA, Ducey MJ, Beers TW, & Husch B. (2016).Forest mensuration. Wiley.
Khattree R, Naik DN (1999) Applied multivariate statistics with SAS software, 2nd edn. SAS Institute Inc., Cary
Krisnawati H, Wang Y, Ades PK (2010) Generalized heightdiameter models for Acacia mangium willd. plantations in South Sumatra. Indonesian J For Res. 7(1):1–19. https://doi.org/10.20886/ijfr.2010.7.1.119
Kutchartt E, Gayoso J, Pirotti F, Bucarey Á, Guerra J, Hernández J, Corvalán P, Drápela K, Olson M, Zwanzig M (2021) Aboveground tree biomass of Araucaria araucana in southern Chile: measurements and multiobjective optimization of biomass models. iForest 14(1):61–70. https://doi.org/10.3832/ifor3492013
Kutchartt E, Hernández J, Corvalán P, Promis Á, Pirotti F (2022) Detecting and evaluating disturbance in temperate rainforest with Sentinel2, machine learning and forest parameters. Int Arch Photogramm Remote Sens Spatial Inf Sci. XLIIIB32022:913–920. https://doi.org/10.5194/isprsarchivesXLIIIB320229132022
Lexerød NL, Eid T (2006) An evaluation of different diameter diversity indices based on criteria related to forest management planning. For Ecol Manage 222(1–3):17–28. https://doi.org/10.1016/j.foreco.2005.10.046
Magnin A, Torres C, Stecconi M, Villalba R, Puntieri J (2021) Influence of trunk forking on height and diameter growth in an evenaged stand of Nothofagus pumilio. NZ J Bot 60:45–59. https://doi.org/10.1080/0028825X.2021.1920433
Marchi M (2019) Nonlinear versus linearised model on stand density model fitting and stand density index calculation: analysis of coefficients estimation via simulation. J For Res 30(5):1595–1602. https://doi.org/10.1007/s11676019009670
Mehtätalo L, deMiguel S, Gregoire TG (2015) Modeling heightdiameter curves for prediction. Canadian Journal of Forest Research 45(7):826–837. https://doi.org/10.1139/cjfr20150054
Mehtätalo L, Kansanen K (2020) lmfor: functions for forest biometrics. R package version 1.5. https://CRAN.Rproject.org/package=lmfor
Molina J, Martín A, Drake F, Martín L, Herrera M (2015) Fragmentation of Araucaria araucana forests in Chile: quantification and correlation with structural variables. iForest 9(2):244–252. https://doi.org/10.3832/ifor1399008
Montaldo P (1974) La bioecologia de Araucaria araucana (Mol.) Koch. Inst. Forestal LatinoAmericano, Bol. Tecn. 46.
Mundo IA, Kitzberger T, Roig Juñent FA, Villalba R, Barrera D (2013) Fire history in the Araucaria araucana forests of Argentina: human and climate influences. Int J Wildland Fire 22:194–206. https://doi.org/10.1071/WF11164
Ogana FN (2021) A mixedeffects heightdiameter model for Gmelina arborea Roxb stands in southwest Nigeria. J for Res 27:1–7. https://doi.org/10.1080/13416979.2021.1989131
Özçelik R, Yavuz H, Karatepe Y, Gürlevik N, Kiriş R (2014) Development of ecoregionbased heightdiameter models for 3 economically important tree species of southern Turkey. Turk J Agric for 38(3):399–412
Özçelik R, Cao QV, Trincado G, Göçer N (2018) Predicting tree height from tree diameter and dominant height using mixedeffects and quantile regression models for two species in Turkey. For Ecol Manage 419420:240–248. https://doi.org/10.1016/j.foreco.2018.03.051
Paulo JA, Tomé J, Tomé M (2011) Nonlinear fixed and random generalized height–diameter models for Portuguese cork oak stands. Ann for Sci 68(2):295–309. https://doi.org/10.1007/s135950110041y
Pinheiro J, Bates D (2000) Mixedeffects models in S and SPLUS. Springer, New York
Pinheiro J, Bates D, DebRoy S, Sarkar D, R Core Team (2021) nlme: linear and nonlinear mixed effects models. R package version 3.1–153. R package version 3.1–153.
Premoli A, Quiroga P, Gardner M (2013) Araucaria araucana. The IUCN Red List of Threatened Species 2013: e.T31355A2805113. https://doi.org/10.2305/IUCN.UK.20131.RLTS.T31355A2805113.en
Pretzsch H (2010) Reevaluation of allometry: stateoftheart and perspective regarding individuals and stands of woody plants. In: Lüttge et al (eds) Progress in Botany 71. pp 339–369. https://doi.org/10.1007/9783642021671_13
Raptis DI, Kazana V, Kazaklis A, Stamatiou C (2021) Mixedeffects height–diameter models for black pine (Pinus nigra Arn.) forest management. Trees 35(4):1167–1183. https://doi.org/10.1007/s0046802102106x
R Core Team (2021) R: a language and environment for statistical computing. R Foundation for Statistical Computing, Vienna. URL https://www.Rproject.org/
Rebertus AJ, Veblen TT (1993) Structure and treefall gap dynamics of oldgrowth Nothofagus forests in Tierra del Fuego. Argentina J Veget Sci 4(5):641–654. https://doi.org/10.2307/3236129
Reineke LH (1933) Perfecting a standdensity index for evenaged forests. J Agric Res 46:627–638
Schmidt M, Kiviste A, von Gadow K (2011) A spatially explicit height–diameter model for Scots pine in Estonia. Eur J Forest Res 130(2):303–315. https://doi.org/10.1007/s1034201004348
Sharma RP, Vacek Z, Vacek S, Kučera M (2019) Modelling individual tree heightdiameter relationships for multilayered and multispecies forests in central Europe. Trees 33(1):103–119. https://doi.org/10.1007/s0046801817624
Sileshi GW (2014) A critical review of forest biomass estimation models, common mistakes and corrective measures. For Ecol Manage 329:237–254. https://doi.org/10.1016/j.foreco.2014.06.026
Subedi MR, Oli BN, Shrestha S, Chhin S (2018) Heightdiameter modeling of Cinnamomum tamala grown in natural forest in midhill of Nepal. Int J For Res 2018, 6583948. https://doi.org/10.1155/2018/6583948
Tarmu T, Laarmann D, Kiviste A (2020) Mean height or dominant height–what to prefer for modelling the site index of Estonian forests? For Stud 72(1):121–138. https://doi.org/10.2478/fsmu20200010
Tischer A, Zwanzig M, Frischbier N (2020) Spatiotemporal statistics: analysis of spatially and temporally correlated throughfall data: exploring and considering dependency and heterogeneity. In: Levia et al (eds) ForestWater Interactions. Ecological Studies. Springer, Cham. https://doi.org/10.1007/9783030260866_8
Trincado G, Van der Schaaf CL, Burkhart HE (2007) Regional mixedeffects height–diameter models for loblolly pine (Pinus taeda L.) plantations. Eur J For Res. 126(2):253–262. https://doi.org/10.1007/s1034200601417
Veblen TT (1982) Regeneration patterns in Araucaria araucana forests in Chile. J Biogeogr 9(1):11–28. https://doi.org/10.2307/2844727
Vonesh E, Chinchilli VM (1997) Linear and nonlinear models for the analysis of repeated measurements. Marcel Dekker Inc, New York. https://doi.org/10.1201/9781482293272
Von Gadow K (2005) Forsteinrichtung: analyse und entwurf der Waldentwicklung. Universitätsverlag Göttingen.
Williams RA (1996) Stand density index for loblolly pine plantations in North Louisiana. South J Appl for 20(2):110–113. https://doi.org/10.1093/sjaf/20.2.110
Wykoff WR (1990) A basal area increment model for individual conifers in the Northern Rocky Mountains. For Sci 36:1077–1104. https://doi.org/10.1093/forestscience/36.4.1077
Xie L, Widagdo RA, Dong L, Li F (2020) Modeling height–diameter relationships for mixedspecies plantations of Fraxinus mandshurica Rupr. and Larix olgensis Henry in northeastern China. Forests. 11(6):610. https://doi.org/10.3390/f11060610
Yuancai L, Parresol BR (2001) Remarks on heightdiameter modeling. Research Note SRS10. US Department of Agriculture, Forest Service, Southeastern Research Station. p. 8.
ZamoranoElgueta C, Cayuela L, GonzalezEspinosa M, Lara A, ParraVazquez MR (2012) Impacts of cattle on the South American temperate forests: challenges for the conservation of the endangered monkey puzzle tree (Araucaria araucana) in Chile. Biol Cons 152:110–118. https://doi.org/10.1016/j.biocon.2012.03.037
Zeide B (1995) A relationship between size of trees and their number. For Ecol Manage 72(2–3):265–272
Zhang L (1997) Crossvalidation of nonlinear growth functions for modelling tree heightdiameter relationships. Ann Bot 79(3):251–257. https://doi.org/10.1006/anbo.1996.0334
Zhou X, Kutchartt E, Hernández J, Corvalán P, Promis Á, Zwanzig M (2022) Tree stem diameter and height of Araucaria araucana, Nothofagus pumilio and Nothofagus dombeyi in mixed stands affected to different levels by anthropogenic disturbance in southcentral Chile. Zenodo. https://doi.org/10.5281/zenodo.7411420
Zhou X, Zwanzig M (2022) Steps for tree height model development and calibration with R. Zenodo. https://doi.org/10.5281/zenodo.7411868
Zu X, Li Q, Ni C, Qin X, Nigh G (2016) Analysis and comparison of combinations among fitting NLME and predictors of random parameters and response variables. Scientia Silvae Sinicae 52(10):72–79. https://doi.org/10.11707/j.10017488.20161009
Zwanzig M, Schlicht R, Frischbier N, Berger U (2020) Primary steps in analyzing data: tasks and tools for a systematic data exploration. In: Levia DF, CarlyleMoses DE, Iida S, Michalzik B, Nanko K, Tischer A (eds) ForestWater Interactions. Ecological Studies vol 240. Springer, Cham. https://doi.org/10.1007/9783030260866_7
Acknowledgements
We thank Valentina González, Daniel Burger, Javiera Aldea, Paula Sandoval, Martin Lotina, Javier Hernández, Ana María Acuña, Bernardita Navarrete, and Joel Hernández for their contribution in the field measurements. Our gratitude to Dr. Javier Guerra for his contribution in the prospection work. A special thanks to the local foresters Ing. Jaime Videla and Ing. Leonardo Araya for the authorization and guidance in the protected areas.
Code availability
The custom Rcode and software application generated during the current study are available as GitHub repository published on Zenodo: https://doi.org/10.5281/zenodo.7411868
Funding
Open Access funding enabled and organized by Projekt DEAL. XZ and MZ were funded by the Federal Ministry of Education and Research (BMBF) and the Freestate of Saxony under the Excellence Strategy of the Federal Government and the Länder. The field data collection was funded through the project 016/2019 “Indicadores fenológicos y estructurales de alteración de hábitat en bosques de Araucaria,” being part of the Fondo de Investigación del Bosque Nativo (FIBN) of the Corporación Nacional Forestal (CONAF) and the Ministry of Agriculture of Chile. The second author is supported by CONICYT doctoral scholarship.
Author information
Authors and Affiliations
Contributions
Conceptualization, JH, PC, and AP conceived the field study, and XZ and MZ designed the methodological approach; methodology, EK, JH, PC, and AP performed the empirical observations and organized the data base; formal analysis and investigation, XZ developed the models, and XZ, EK, and MZ interpreted the analytical results; writing — original draft preparation, XZ prepared the original draft; writing — review and editing, MZ reviewed and revised the original draft; funding acquisition, JH, PC, AP, and MZ acquired funds; supervision, MZ and EK supervised XZ. The authors read and approved the final manuscript.
Corresponding author
Ethics declarations
Ethics approval and consent to participate
The authors declare that they obtained the permission of the Department of Protected Wildlife Areas — CONAF to get access to the National Park Conguillio and National Reserve Malalcahuello by Authorization No. 05/2019 IX.
Consent for publication
All authors gave their informed consent to this publication and its content.
Competing interests
The authors declare that they have no competing interests.
Additional information
Handling editor: Shuguang Liu
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Appendices
Appendices
1.1 Detailed workflow for heightdiameter model development
1.1.1 General principle of heightdiameter model development
The definition of appropriate heightdiameter models can be a challenging task. Here, we have taken up recommendations by Tischer et al. (2020). We begin with an exploration of simple models and continue by examining model extensions that account for covariates and random effects. An R script for this workflow is presented by Zhou & Zwanzig (2022).
1.1.2 Simple models selection
The selection of H–D models requires consideration of several characteristics such as: desirable mathematical properties (e.g. number of parameters), possible biological interpretations of the parameters and satisfactory predictions of the heightdiameter relationships (Yuancai and Parresol 2001; Krisnawati et al. 2010). To determine the most appropriate simple model for A. araucana and N. pumilio, 16 nonlinear models that have been used in many studies (Mehtätalo et al. 2015) were chosen as candidates (Table 12).
They can all be written in the following general form (Adame et al. 2008):
where \({H}_{i}\) is the \(i\)th observation of tree height (m), \(DB{H}_{i}\) is the \(i\)th observation of diameter at breast height (cm), \(\phi\) is the vector of parameters to be estimated, \({\varepsilon }_{i}\) is the random error term and \(i\) is the \(i\)th observation with \(i=\mathrm{1,2},\dots ,n\).
Among the 16 candidate models are seven nonlinear models with two parameters and nine nonlinear models with three parameters. Nonlinear models containing four parameters were not involved in this study to prevent overparameterization and overcomplication of the models. The constant height at breast height (1.3) is included in the right side of all H–D models, thus avoiding negative height estimates for small trees (Krisnawati et al. 2010).
The performance of the different models was assessed with four types of criteria: (1) evaluation of uncertainty and identifiability of the parameter estimates with percentage relative standard error (PRSE); (2) assessment of the goodnessoffit by the root mean square error (RMSE), mean absolute error (MAE) and the coefficient of determination (R^{2}); (3) evaluation of the prediction ability of these simple models on fitting database based on the mean RMSE of the tenfold CrossValidation repeated 10 times (Ciceu et al. 2020), and (4) visual analysis based on studentized residuals plots (Huang et al. 1992).
The expressions of these statistics are summarized as follows:
where \(P\) is the parameter estimate, \(SE\left(P\right)\) is the parameter standard error, \({H}_{i}\) and \({\widehat{H}}_{i}\) are the observed and predicted values of the height of tree \(i\), \(\overline{H}\) is the average value of the height and \(n\) is the number of trees.
In general, lower RMSE and MAE values and higher R^{2} values indicated that the model was providing a better goodnessoffit, whereas lower PRSE values indicated a proper identifiability of the parameter estimates. In many practical applications, parameter estimates were considered unreliable when PRSE exceed 25 30%—a rule of thumb for PRSE reported by Sileshi (2014). For each alternative model, the maximum PRSE value of its parameters would be taken as the PRSE value for that model and when this value exceeded 25%, the model would be rejected.
The nonlinear H–D models were fitted using the nlsfunction of the Rpackage ‘stats’ (R Core Team 2021). A set of the initial values of the parameters were obtained using the startHDmodelfunction of the Rpackage ‘lmfor’ (Mehtätalo and Kansanen 2020) in order to avoid inappropriate initial value settings from causing the model to fail to fit on the fitting database. To ensure that each fitted simple model for given height and DBH data was globally optimal rather than locally optimal, multiple different sets of initial values were also assumed and tested and the Gauss–Newton algorithm that comes with the nlsfunction was used to determine the nonlinear least squares estimates of the parameters.
The assumption of nonlinear least squares is that the error terms are independent and identically distributed with zero mean and constant variance (Huang et al. 1992). Therefore, after determining the best simple model for each species, if the residual diagnostics showed heteroskedasticity, it would be improved by weighted nonlinear least squares (WLS). According to Huang et al. (1992), the weighting factor could be set as:
where the alternative values of \(k\) are 0.5, 1, 1.5, 2, 2.5, 3. After the simple model with the best fit was identified, these weighting factors were tested for further model improvement.
1.1.3 Stand variables
In order to study the variability of the relationship between tree height and diameter at stand level, a series of variables based on the original data available were selected as alternative covariates to be added to the simple model in subsequent step. These 25 variables were classified into six main categories, namely tree, stand, diameter diversity, density and competition, geographic information, and mixed stand (see Table 17 in Appendix).
The stand variables describing the tree level were mean DBH (DBH_{m}) and mean tree height (H_{m}).
The quadratic mean DBH and tree height (Dq and Hq), maximum DBH and tree height (DBH_{max} and H_{max}), range of DBHs and tree heights (range DBH and range H), height of the tree with the largest DBH (H_{Dmax}), range of heights of the trees with the smallest and the largest DBH (range H_{D}), and dominant diameter and tree height (DBH_{dom} and H_{dom}) were used to describe stand level. H_{dom} was reduced to the mean height of the three thickest trees in the stand, although a more common definition in forestry research is to consider dominant height as the average height of the 100 thickest trees per hectare (see Tarmu et al. 2020). The simplification made here was on the one hand because only three trees in some stands in the database were measured for height and on the other hand the consideration that the inclusion of covariates requiring many premeasured tree heights in the model for predicting tree height would significantly reduce the applicability of the model. For the reason of standardization, the calculation of the dominant diameter has been simplified in the same way.
Diameter diversity was described with the help of the Shannon index (Hd), the Gini coefficient (GCd), and the coefficient of variation of DBH (Vard). Where the Shannon index was originally used to describe species diversity, it was later adapted as an index of diameter diversity by replacing the number of species with the number of diameter classes of individual trees (Dănescu et al. 2016). The Gini coefficient was considered because it not only did not require subjective assignment of classes of diameters but also measured heterogeneity (Lexerød and Eid 2006). To calculate the Gini coefficient, all trees within the sampling unit were arranged in ascending order of diameter, and the result was between zero and one, where the smallest value of zero represented the fact that all trees were the same size. The coefficient of variation of DBH could help to compare the dispersion of diameters in stands with different mean diameters.
Tree growth was often exceptionally sensitive to the differences of stand density because the density would affect the spatial distribution of the climate within a stand, especially light and temperature. The most common measures of stand density include, but are not limited to, the number of trees per hectare, basal area per hectare (BA), and stand density index, each of which has its advantages and disadvantages, and there is no universally applicable measure of stand density (Von Gadow 2005). For example, the number of trees per hectare (N), although very simple and straightforward, lacked the second dimension of individual size, and it was not suitable for describing stand density anymore if the mean size of the population was different (Zeide 1995). The stand density index (SDI) was introduced by Reineke in 1933 (Reineke 1933) and is now widely used in the USA (e.g., Marchi 2019). The theoretical SDI is an abstract amount indicating the expected number of trees when Dq is 25 cm. When using SDI, the Reineke constant (1.605) is often used directly, but its generality has actually been questioned (e.g., Williams 1996). The basal area of larger trees (BAL), which originally referred to the sum of the basal areas of trees larger than the diameter of a given tree (Wykoff 1990), was also simplified in this study to the sum of the basal areas of the three largest diameter trees in the stand. The simplification was since there was not a large amount of data available to determine the most reasonable reference diameter. BAL was an important explanatory variable in many modeling approaches (e.g., Xie et al. 2020), but species blindness needed to be taken into account when using BAL indices in mixed stands (Von Gadow 2005).
Geographic information included altitude (ALT), latitude, and longitude (lat and lng).
Mixed stand was simply considered as the proportion of the basal area of each species in relation to the total basal area (RBA) and the proportion of the number of each species in relation to the total number of trees (RN), as the original database did not include information on the distance between trees.
Some of the stand variables considered both species involved and species uninvolved.
To see if specific quantitative standlevel properties are characteristic for the four levels of anthropogenic disturbance, a kmeans cluster analysis with four predefined clusters was conducted on the stand variables involved in this study. For this, the “kmeans” function of Rpackage “stats” (R Core Team 2021) was used. The visualization of the clustering results used the “fviz_cluster”—function of Rpackage “factoextra” (Kassambara and Mundt 2020), which represents observations in a plot based on principal components. The detailed results of the principal components analysis were obtained using prcomp() of the “stats” package (R Core Team 2021).
1.1.4 Inclusion of stand variables in the models
Simple H–D models could be extended upwards to generalized H–D models by including stand variables, thus expressing the variability detected between stands (Adame et al. 2008; Mehtätalo et al. 2015) and improving model fitting and prediction (Sharma et al. 2019). Three different approaches of including stand variables were tested. The first two represented a twostage approach, which was introduced by Ferguson and Leech (1978) and was still often applied (e.g. Ciceu et al. 2020). In the first stage, separate heightdiameter curves were fitted for each stand for both approaches. In the second stage, the first approach was to progressively add to the simple model the stand variables with the highest correlation to the model parameters. The second approach was using best subset regression to search for the best linear fit for different numbers of stand variables for each parameter. The third “exhaustive” approach was to explore all possible combinations to find the stand variable or combination of stand variables that gave the largest improvement to the goodnessoffit statistics, while limiting the total number of stand variables to be included.
In the first approach, the Pearson correlation coefficient was used. Therefore, the calculated stand variables would first be checked for the normal distribution, and a simple arithmetic transformation (log, inverse, and squared, Sharma et al. 2019) would be carried out for the stand variables that did not fit a normal distribution. If the transformed stand variables still did not conform to a normal distribution, the Spearman correlation coefficient was used as an alternative. RD was not considered in the first approach because each tree had a different RD, while the other stand variables were the same for all the trees in the same stand and the inclusion of RD would overincrease the significance of the correlation. The rule for interpreting the correlation coefficient referred to Hinkle et al. (2003).
It is important to note that because one of the purposes of developing the heightdiameter model is to predict tree height, stand variables that required large amounts of tree height data for calculation were excluded (e.g., H_{m}). Since most of the stand variables considered were calculated based on DBH, the problem of collinearity might arise when more than one stand variable is included. This problem was assessed by calculating variance inflation factor (VIF) and considering that, as a rule of thumb, a VIF value that exceeds 5 or 10 indicates a problematic amount of collinearity (Gareth et al. 2013).
1.1.5 Comparison of H–D models among different sites
Because of the apparent multiple layers of nested relationships (sitestandplot) in the database, this spatial correlation in the data must be considered when developing the H–D models. To compare the differences between the heightdiameter relationships at the four sites of our study area, the nonlinear extra sum of squares Ftest (Bates and Watts 1988; Huang et al. 2000) and the test introduced by Lakkis and Jones (Khattree and Naik 1999) were used. For this purpose, both full and reduced models were required. The full model had a different set of parameters for each site, while the reduced model shared the same set of parameters for all four sites. In this study, the reduced model is the H–D model developed in the previous step including stand variables. The full model would be extended on this by the indicator variable approach (Bates and Watts 1988), i.e., extending each parameter to include an associated parameter as well as a dummy variable (e.g., \({\Delta }_{1}{x}_{1}\) for La Fusta) to distinguish among sites (Adame et al. 2008).
Since there are 4 different sites, 3 dummy variables (\({x}_{1}\) to \({x}_{3}\)) were needed, which were defined as follows:

If \(site=LaFusta,{x}_{1}=1,{x}_{2}=0,{x}_{3}=0\)

If \(site=Conguillio,{x}_{1}=0,{x}_{2}=1,{x}_{3}=0\)

If \(site=Malalcahuello,{x}_{1}=0,{x}_{2}=0,{x}_{3}=1\)

If \(site=ElNaranjo,{x}_{1}=0,{x}_{2}=0,{x}_{3}=0\)
The equality of full model and reduced model was tested by considering the null hypothesis that all associated parameters were equal and equal to zero (Huang et al. 2000).
The expressions of the Ftest and LakkisJones test mentioned above are as follows:
where \(SS{E}_{R}\) and \(SS{E}_{F}\) are the error sum of squares of the reduced model and full model, \(d{f}_{R}\) and \(d{f}_{F}\) are the degrees of freedom of the reduced and full model, and \(n\) is the number of observations.
The null hypothesis was accepted only if both tests passed (\(F\le {F}_{critical}\left(1\alpha ,d{f}_{R}d{f}_{F},d{f}_{F}\right),\alpha =0.05\) and \(L\le {L}_{critical}\left(1\alpha ,d{f}_{R}d{f}_{F}\right),\alpha =0.05\)), i.e., the reduced model was applied to all four sites.
Both the reduced model and the full model were fitted with Rfunction “nls.” The residuals are needed to be diagnosed again to check for heteroskedasticity.
1.1.6 Nonlinear mixedeffects models
The effect of different stands on the heightdiameter relationship was addressed using nonlinear mixedeffects models. The general form of the nonlinear mixedeffects H–D model can be written as follows (Pinheiro and Bates 2000):
where \(M\) is the number of stands, \({n}_{i}\) is the number of trees in the \(i\) th stand, \(f\) is a general, realvalued, differentiable function of a standspecific parameter vector \({\phi }_{ij}\), \(DB{H}_{ij}\) is a covariate vector, and \({\varepsilon }_{ij}\) is a normally distributed withinstand error term.
An important characteristic of the nonlinear mixedeffects models was that their parameters were divided into two groups: fixed effects parameters, which showed trends in tree heights common to all stands, and random effects parameters, which represented differences among stands. The parameter vector \({\phi }_{ij}\) can therefore be defined as follows (Pinheiro and Bates 2006):
where \(\beta\) is a pdimensional vector of fixed effects, \({b}_{i}\) is a qdimensional random effects vector associated with the \(i\) th stand with variance–covariance matrix \(D\), and \({A}_{ij}\) and \({B}_{ij}\) are the design matrix for fixed and random effects.
Taking the example of Fang and Bailey (2001) and Dorado et al. (2006), the development of the mixedeffects models in this study involved three main steps.
The first step was to determine the parameter effects, which was to specify the nature of the parameters of the model as fixed and random effects or purely fixed effects. Pinheiro and Bates (2006) suggested that, when convergence was feasible, all parameters in the model were first considered to including both random and fixed effects. In this study, various variants of the mixedeffects models formed by all possible combinations of random and fixed effects parameters were fitted to the fitting database. When different models were fitted to the same database, comparisons between models could be either nested models or nonnested models. Nested models were compared using a likelihood ratio test (LRT). For nonnested models, the LRT test was not available, in which case the models were compared using the Akaike information criterion (AIC). A smaller value of AIC indicated a better fit of the model to the database.
The second step was to determine the withinstand variance–covariance structure to account for the variability between tree heights in the same stand. To achieve this, it was important to address both the heteroskedasticity and autocorrelation structure components (Fang and Bailey 2001). Residual diagnostics were used to determine if there was heteroskedasticity, and if so, one of the three common variance functions (power function, exponential function, and constant plus power function) with the greatest reduction in heteroskedasticity would be used to resolve this problem. As there were no repeated measurements in the database, it was reasonable to assume that there was no temporal correlation between these observations. Although trees originating from the same stand all shared the same environment, trees were randomly selected when measuring height, so it was also reasonable to assume that these observations were not clumped in the same stand, i.e., there was no pattern of spatial correlation.
Dorado et al. (2006) considered that the expression for the withinstand variance–covariance matrix has the following special structure, when no patterns of temporal or spatial correlation among observations are shown in the data and only a weighting factor for balancing the error variance is considered:
where \({\sigma }^{2}\) is a scaling factor for the error dispersion given by the value of the residual variance of the model and \({G}_{i}\) is the diagonal matrix describing the nonconstant variance.
The third step was to determine the structure of the amongstand variance–covariance matrix \(D\), which was common to all stands and defined the variability that existed among stands. A general \(D\) with \(n\) parameters has the following equation:
where \({\sigma }_{1}^{2}\) and \({\sigma }_{n}^{2}\) are the variance of the first and the \(n\) th parameters and \({\sigma }_{1}{\sigma }_{n}\) and \({\sigma }_{n}{\sigma }_{1}\) are the mutual covariances of the first and the \(n\) th parameters.
The nonlinear mixedeffects H–D models were fitted using function nlme() of package “nlme” (Pinheiro et al. 2021). The algorithm used was loglikelihood maximization (ML).
1.1.7 Calibration of the nonlinear mixedeffects models and height prediction
Using the nonlinear mixedeffects models developed, tree height can be predicted in stands to which the model has not been previously fitted. There are two types of prediction, a fixed effects response pattern, where only diameters and the stand variables involved in the model are measured in the new stand and a calibrated response pattern, where not only diameters and stand variables but also the height of a subsample of trees in the new stand are measured (Fang and Bailey 2001).
In the first case, the expected value 0 was used for all random parameters, and the prediction of the fixed part was the standardized generalized heightdiameter curve (Dorado et al. 2006). In the second case, the subsample of tree heights can be used to predict the random effects vector \({b}_{i}\) for the new stand, using the empirical best unbiased predictor (EBLUP) technique (Vonesh and Chinchilli 1997), expressed as follows:
where \(\widehat{D}\) is the amongstand variance–covariance matrix, common for all stands, and estimated in the general fitting of the model, \({\widehat{R}}_{i}\) is the withinstand variance–covariance matrix, \({\widehat{Z}}_{i}\) is the vector of partial derivatives with respect to \({b}_{i}\), and \({\widehat{\varepsilon }}_{i}\) is the residual vector given by the difference between the observed height value and the value predicted using the model including only fixed effects. The superscript \(T\) is the matrix transpose operator, and \(1\) is the inverse matrix.
For this calibrated response pattern, this study evaluated alternatives with different sampling sizes within each stand, i.e., a random selection of one to eight sample trees for height. Height for other trees in the stand can be predicted with the help of the following equation after determining the random effects in the new stand:
It is important to note that EBLUP was involved in the process of linearizing a nonlinear model using a firstorder Taylor approximation, and the linearized base point could be either 0 or the iterative final value. In forestry, 0 is generally applied as the linearized base point to predict random effects (Eq. 16), so it was preferable to use the same base point for linear prediction of the tree height (Eq. 17). This was because the accuracy of the prediction is substantially reduced when the base point of predicted random effect differs from the prediction of tree height (Zu et al. 2016).
To assess the prediction performance of the developed nonlinear mixedeffects models, data included in validating database were used. The calibrated H–D models were applied to the remaining trees in the same stand for which height measurements were available. Eight calibration design alternatives were evaluated and compared with the predictions obtained from ONLS techniques in the same stand with selected simple and generalized nonlinear H–D models and with a typical fixed effects response pattern.
The prediction performance of the models was evaluated using the following statistical measures (Zu et al. 2016; Sharma et al. 2019) as follows:
where \({H}_{ij}\) and \({\widehat{H}}_{ij}\) are observed and predicted values of the height of the \(j\) th tree in the \(i\) th stand, \({n}_{i}\) is the number of predicted values for the \(i\) th stand, \(m\) is the number of stands involved in the height prediction analysis, and \({\overline{H}}_{i}\) is the mean value of measured heights in the \(i\) th stand.
Random selection of the sample trees was repeated 1000 times to obtain an average estimate for the final results.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Zhou, X., Kutchartt, E., Hernández, J. et al. Determination of optimal tree height models and calibration designs for Araucaria araucana and Nothofagus pumilio in mixed stands affected to different levels by anthropogenic disturbance in SouthCentral Chile. Annals of Forest Science 80, 18 (2023). https://doi.org/10.1186/s13595023011859
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s13595023011859