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 South-Central Chile

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 Araucaria-Nothofagus 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 mixed-effects models.

Context

Araucaria-Nothofagus 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 stand-level properties were calculated to determine forest structure, which was compared among stands by cluster analysis. To identify the optimal height-diameter (H–D) model, simple models including diameter only as well as generalized models including stand variables were tested, each additionally extended by a nonlinear mixed-effects (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 height-diameter relationships, the same H–D model can be used for stands representing different anthropogenic disturbance levels and recalibrated by cost-effective 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 South-Central 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. araucana-N. pumilio forests are mainly found between 1000 and 1600 m asl in the Andes Cordillera, but also mixed A. araucana-N. antarctica forests could be found at the same altitude. Around 48% of the A. araucana dominated and natural forests in Chile are protected as state-protected 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 human-caused 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; Zamorano-Elgueta et al. 2012; Molina et al. 2015; Hernández et al. 2022). Thus, a large proportion of Araucaria-Nothofagus 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; Fuentes-Ramí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 Araucaria-Nothofagus 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 slow-growing relict conifer and a relatively short-lived shade-intolerant hardwood, are both important species with high ecological and economic importance. Technical reports by Gayoso (2013a, 2013b) presented height-diameter models for these two species, but they have not been assessed by cross-validation. 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 stand-specific 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 height-diameter 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 time-consuming, 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 height-diameter 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 height-diameter 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 height-diameter 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 cross-validation-based root-mean-square 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 height-diameter relationships among stands. However, height-diameter relationships vary between stands and are influenced by the growing environment and stand conditions. Even within the same stand, height-diameter relationships can change over time (Schmidt et al. 2011). Generalized H–D models use stand variables to estimate local variation of the height-diameter 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 two-stage approach (Ferguson and Leech 1978) has often been used in later studies (e.g., Ciceu et al. 2020). The two-stage approach involves fitting a separate height-diameter 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 height-diameter 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 mixed-effects (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 height-diameter relationship (Calama and Montero 2004). And NLME technique can be applied both to simple and generalized H–D models (Gómez-Garcí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 South-Central 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 South-Central 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 Araucaria-Nothofagus 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 sandy-loamy 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 cation-exchange capacity is very low and the base saturation only 2–5%.

Fig. 1

Study area and distribution of 12 selected stands in the northeastern sector of the Araucanía region in southern Chile. Stands are located in La Fusta, Conguillío (National Park), Malalcahuello (National Reserve) and El Naranjo and are distinguished by four levels of anthropogenic disturbance: none, low, medium, high

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 long-lived, slow-growing but huge relict conifer native to South-Central 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. araucana-Nothofagus forests (Veblen 1982).

Nothofagus pumilio (Poepp. & Endl.) Krasser is a broad-leaved 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 post-fire (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² per hectare. DBH was measured for all selected trees, and height was measured for one-third 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 height-diameter data pairs of A. araucana and 380 data pairs of N. pumilio were extracted for the study of the height-diameter relationship for these two species. In order to develop a reliable height-diameter 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).

Table 1 Summary of the complete database

Fig. 2

Height plotted against DBH for the complete, fitting, and validating database of A. araucana (a, b, and c) and N. pumilio (d, e, and f)

2.3 Height-diameter 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 best-performing 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 mixed-effects 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.

Fig. 3

Workflow for height-diameter model development and calibration. A detailed description of the approaches and techniques applied during each step of the workflow is provided in the appendix. See Table 12 for a list of simple models (SM1–SM16) tested

A description of the specific technical processes for each step can be found in the annex: Detailed workflow for height-diameter 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).

Table 2 Overview of common and more advanced tree height models developed

3 Results

In the following, first, the results of the k-means 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 k-means 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.

Fig. 4

Clusters of anthropogenic disturbances levels. Results refer to a k-means clustering analysis with observations represented along the first two principal components explaining ca. 33% and 20% of the observed variance in forest structure. Numbers refer to stands with varying levels of anthropogenic disturbance: none (stands 4, 5, and 6), low (stand s7, 8, and 9), medium (stands 10, 11, and 12), and high (stands 1, 2, and 3)

3.2 Height-diameter 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 cross-validation data, the lowest MAE and the highest R² values (Table 3).

Table 3 A. araucana simple model selection results. For formulae of models SM1-SM16, see Appendix Table 12. Prediction statistic refers to the cross-validation (CV)-based prediction accuracy

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.

Table 4 Results from the inclusion of stand variables for SM12 AA

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:

$${H}_{ij}=1.3+\left(a0+a1\times {H}_{dom}A{A}_{i}\right)\times {e}^{-\left(b0+b1\times RangeDBHA{A}_{i}\right)\times {e}^{-c0\times {DBH}_{ij}}}+{\varepsilon }_{ij}$$

(1)

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 height-diameter relationship were simulated (Fig. 5). The effect of ${H}_{dom}AA$ on the height-diameter 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.

Fig. 5

Effects of dominant height of A. araucana (${H}_{dom}AA$) and the difference between the diameters at breast height of A. araucana ($RangeDBHAA$) on the height-diameter relationships of A. araucana. The curves were produced from the parameter estimates obtained when fitting Eq. 1. The values of the stand variables were replaced using the mean value except for the stand variable of interest, which slowly varied from the minimum to the maximum at the same spacing

Comparison of H–D models among different sites

Equation (1) was used as a reduced model when comparing height-diameter relationships among four different sites. The full model extended by the indicator variable method could be written as follows:

$${H}_{ij}=1.3+\left(a0+{\Delta }_{1}{x}_{1}+{\Delta }_{2}{x}_{2}+{\Delta }_{3}{x}_{3}+a1\times {H}_{dom}A{A}_{i}\right)\times {e}^{-\left(b0+{\phi }_{1}{x}_{1}+{\phi }_{2}{x}_{2}+{\phi }_{3}{x}_{3}+b1\times RangeDBHA{A}_{i}\right)\times {e}^{-\left(c0+{\psi }_{1}{x}_{1}+{\psi }_{2}{x}_{2}+{\psi }_{3}{x}_{3}\right)\times {DBH}_{ij}}}+{\varepsilon }_{ij}$$

(2)

Both the F-test and the Lakkis-Jones 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 height-diameter 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 height-diameter 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).

Table 5 Results of statistical tests comparing the height-diameter relationship of A. araucana and N. pumilio among different sites

Nonlinear mixed effects models

First, a nonlinear mixed-effects 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 within-stand variance (${\sigma }^{2}$) did not drop significantly, suggesting the existence of a pattern of variability that cannot be explained by differences among stands.

Table 6 Estimated parameters and goodness of fit of the H–D models for A. araucana, developed on the basis of the fitting database. For estimates for the models fitted to the complete database, see Table 15 in Appendix

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.

Table 7 Prediction bias (%) of the NLME H–D models of A. araucana in stand 2 and stand 12, respectively

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 fixed-effects response pattern and the calibrated response pattern were close, but both were different from the curves obtained from the ONLS.

Fig. 6

Example of the calibration of simple NLME model with five sample trees and generalized NLME model with one sample tree in stands 2 and 12 compared to the fixed effects and ONLS prediction with simple H–D model and generalized H–D model of A. araucana. The ONLS fits of generalized H–D model did not converge in either stand. Therefore, the generalized NLME model was likewise compared with the ONLS prediction of the simple H–D model

3.2.2 Nothofagus pumilio (Poepp. & Endl.) Krasser

Simple model selection

The two-parameter SM1 (Näslund) with the lowest MAE, RMSE, and RMSE (CV) and the highest R² 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.

Table 8 N. pumilio simple model selection results. For formulae of models SM1-SM16, see Table 12. Prediction statistic refers to the cross-validation (CV)-based prediction accuracy

Studentized residuals can be used to identify outliers, which are defined as points that are far from other observations in one-, two- or n-dimensional 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.

Table 9 Results from the inclusion of stand variables for SM1 NP

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:

$${H}_{ij}=1.3+\frac{DB{H}_{ij}^{2}}{{\left(a0+a1\times DBHma{x}_{i}+\left(b1\times D{q}_{i}+b2\times NlogN{P}_{i}\right)\times {DBH}_{ij}\right)}^{2}}+{\varepsilon }_{ij}$$

(3)

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 log-transformed forms of trees per hectare of N. pumilio in the $i$ th stand.

The effects of these three stand variables on the height-diameter relationship were simulated (Fig. 7). The effect of $Dq$ on the height-diameter 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.

Fig. 7

Effects of maximum DBH ($DB{H}_{max}$), quadratic mean DBH ($Dq$), and number of N. pumilio per ha ($NNP$) on the height-diameter relationships of N. pumilio. The curves were produced from the parameter estimates obtained when fitting Eq. 3. The values of the stand variables were replaced using the mean value except for the stand variable of interest, which slowly varied from the minimum to the maximum at the same spacing

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 F-test and the Lakkis-Jones 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 mixed-effects 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 within-stand variance (${\sigma }^{2}$) was even getting higher. This suggested that the variability in the height-diameter 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.

Table 10 Estimated parameters and goodness of fit of the H–D models for N. pumilio, developed on the basis of the fitting database. For estimates for the models fitted to the complete database, see Table 16 in Appendix

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.

Table 11 Prediction bias (%) of the simple NLME H–D model of N. pumilio in stand 2 and stand 12, respectively

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.

Fig. 8

Example of the calibration of simple NLME model with five sample trees and generalized NLME model with one sample tree in stands 2 and 12 compared to the fixed effects and ONLS prediction with simple H–D model and generalized H–D model of N. pumilio. The ONLS fits of generalized H–D model did not converge in either stand. Therefore, the generalized NLME model was likewise compared with the ONLS prediction of the simple H–D 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 South-Central 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 South-Central Chile. As a result of the presented workflow for the development of height-diameter 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 South-Central 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 large-scale environmental factors, model predictions for dependent tree-level 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 two-parameter models and nine three-parameter models. Although convergence is considered one of the challenges of fitting 3-parameter H–D models (Mehtätalo et al. 2015; Ogana 2021), the nine alternative 3-parameter 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 two-parameter models did not suffer from this problem. Except for SM5 (Power), their goodness of fit was in fact very similar to that of the three-parameter models.

The simple H–D model chosen for A. araucana was SM12 (Gompertz), a model whose suitability in describing the height-diameter 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 cross-validation-based RMSE for the accuracy in predicting new data facilitates the identification of accurate and robust models. Many of the height-diameter 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 height-diameter 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ómez-Garcí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 x-axis. 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 even-aged or even-sized 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 counter-intuitive, 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 site-based 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 mixed-effects 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 fixed-effects response pattern and a calibrated response pattern. The fixed-effects 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 fixed-effects response pattern of both simple and generalized NLME models in even-aged 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 among-stand variance and higher within-stand 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 mixed-effects 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 height-diameter 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 Araucaria-Nothofagus forests.

Appendices

Appendices

1.1 Detailed workflow for height-diameter model development

1.1.1 General principle of height-diameter 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. 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 height-diameter 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).

Table 12 List of 16 simple models tested. Models were classified into 2-parameter and 3-parameter functions, following Mehtätalo et al. (2015), which provides a full list of references for the model origins

They can all be written in the following general form (Adame et al. 2008):

$${H}_{i}=f\left(DB{H}_{i},\phi \right)+{\varepsilon }_{i}$$

(4)

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 over-parameterization and over-complication 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 goodness-of-fit by the root mean square error (RMSE), mean absolute error (MAE) and the coefficient of determination (R²); (3) evaluation of the prediction ability of these simple models on fitting database based on the mean RMSE of the tenfold Cross-Validation 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:

$$PRSE=100\times \frac{SE\left(P\right)}{\left|P\right|}$$

(5)

$$RMSE=\sqrt{\frac{1}{n}{\sum }_{i=1}^{n}{\left({H}_{i}-{\widehat{H}}_{i}\right)}^{2}}$$

(6)

$$MAE={\sum }_{i=1}^{n}\frac{\left|{H}_{i}-{\widehat{H}}_{i}\right|}{n}$$

(7)

$${R}^{2}=1-\frac{{\sum }_{i=1}^{n}{\left({H}_{i}-{\widehat{H}}_{i}\right)}^{2}}{{\sum }_{i=1}^{n}{\left({H}_{i}-\overline{H}\right)}^{2}}$$

(8)

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² values indicated that the model was providing a better goodness-of-fit, 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 nls-function of the R-package ‘stats’ (R Core Team 2021). A set of the initial values of the parameters were obtained using the startHDmodel-function of the R-package ‘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 nls-function 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:

$${w}_{i}=\frac{1}{DB{H}_{i}^{k}}$$

(9)

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 pre-measured 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 stand-level properties are characteristic for the four levels of anthropogenic disturbance, a k-means cluster analysis with four predefined clusters was conducted on the stand variables involved in this study. For this, the “kmeans” function of R-package “stats” (R Core Team 2021) was used. The visualization of the clustering results used the “fviz_cluster”—function of R-package “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 two-stage 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 height-diameter 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 goodness-of-fit 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 over-increase 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 height-diameter 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 (site-stand-plot) in the database, this spatial correlation in the data must be considered when developing the H–D models. To compare the differences between the height-diameter relationships at the four sites of our study area, the nonlinear extra sum of squares F-test (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 F-test and Lakkis-Jones test mentioned above are as follows:

$$F=\frac{\frac{SS{E}_{R}-SS{E}_{F}}{d{f}_{R}-d{f}_{F}}}{\frac{SS{E}_{F}}{d{f}_{F}}}$$

(10)

$$L=2ln\left({\left(\frac{SS{E}_{R}}{SS{E}_{F}}\right)}^\frac{n}{2}\right)$$

(11)

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 R-function “nls.” The residuals are needed to be diagnosed again to check for heteroskedasticity.

1.1.6 Nonlinear mixed-effects models

The effect of different stands on the height-diameter relationship was addressed using nonlinear mixed-effects models. The general form of the nonlinear mixed-effects H–D model can be written as follows (Pinheiro and Bates 2000):

$${H}_{ij}=f\left({\phi }_{ij},DB{H}_{ij}\right)+{\varepsilon }_{ij};i=1,\dots ,M;j=1,\dots ,{n}_{i}$$

(12)

where $M$ is the number of stands, ${n}_{i}$ is the number of trees in the $i$ th stand, $f$ is a general, real-valued, differentiable function of a stand-specific parameter vector ${\phi }_{ij}$, $DB{H}_{ij}$ is a covariate vector, and ${\varepsilon }_{ij}$ is a normally distributed within-stand error term.

An important characteristic of the nonlinear mixed-effects 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):

$${\phi }_{ij}={A}_{ij}\beta +{B}_{ij}{b}_{i};{b}_{i}N\left(0,D\right)$$

(13)

where $\beta$ is a p-dimensional vector of fixed effects, ${b}_{i}$ is a q-dimensional 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 mixed-effects 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 mixed-effects 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 non-nested models. Nested models were compared using a likelihood ratio test (LRT). For non-nested 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 within-stand 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 within-stand 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:

$${R}_{i}={\sigma }^{2}{G}_{i}$$

(14)

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 among-stand 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:

$$D=\left[\begin{array}{ccc}{\sigma }_{1}^{2}& \cdots & {\sigma }_{n}{\sigma }_{1}\\ \vdots & \ddots & \vdots \\ {\sigma }_{1}{\sigma }_{n}& \cdots & {\sigma }_{n}^{2}\end{array}\right]$$

(15)

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 mixed-effects H–D models were fitted using function nlme() of package “nlme” (Pinheiro et al. 2021). The algorithm used was log-likelihood maximization (ML).

1.1.7 Calibration of the nonlinear mixed-effects models and height prediction

Using the nonlinear mixed-effects 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 height-diameter 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:

$${\widehat{b}}_{i}=\widehat{D}{\widehat{Z}}_{i}^{T}{\left({\widehat{R}}_{i}+{\widehat{Z}}_{i}\widehat{D}{\widehat{Z}}_{i}^{T}\right)}^{-1}{\widehat{\varepsilon }}_{i}$$

(16)

where $\widehat{D}$ is the among-stand variance–covariance matrix, common for all stands, and estimated in the general fitting of the model, ${\widehat{R}}_{i}$ is the within-stand 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:

$${H}_{i}\approx f\left({A}_{i}\widehat{\beta },{DBH}_{i}\right)+\widehat{{Z}_{i}}\times \widehat{{b}_{i}}$$

(17)

It is important to note that EBLUP was involved in the process of linearizing a nonlinear model using a first-order 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 mixed-effects 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:

$$MSPE=\frac{1}{{\sum }_{n=1}^{m}{n}_{i}}{\sum }_{i=1}^{m}{\sum }_{j=1}^{{n}_{i}}{\left({H}_{ij}-{\widehat{H}}_{ij}\right)}^{2}$$

(18)

$$MPE=\frac{100}{{\sum }_{n=1}^{m}{n}_{i}}{\sum }_{i=1}^{m}{\sum }_{j=1}^{{n}_{i}}\frac{{H}_{ij}-{\widehat{H}}_{ij}}{{H}_{ij}}$$

(19)

$$MAPE=\frac{1}{{\sum }_{n=1}^{m}{n}_{i}}{\sum }_{i=1}^{m}{\sum }_{j=1}^{{n}_{i}}\frac{\left|{H}_{ij}-{\widehat{H}}_{ij}\right|}{{H}_{ij}}$$

(20)

$$bias\left(\%\right)=\frac{100{\overline{e}}_{i}}{{\overline{H}}_{i}}with{\overline{e}}_{i}={\sum }_{j=1}^{{n}_{i}}\frac{\left({H}_{ij}-{\widehat{H}}_{ij}\right)}{{n}_{i}}$$

(21)

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.

Table 13 Calibration results of NLME H–D models of A. araucana

Table 14 Calibration results of NLME H–D models of N. pumilio. Bold: statistics of selected models

Table 15 Estimated parameters and goodness of fit of the H–D models for A. araucana fitted to the complete database

Table 16 Estimated parameters and goodness of fit of the H–D models for N. pumilio fitted to the complete database

Table 17 Summary table of all stand variables involved (modified based on Ciceu et al. 2020; Dănescu et al. 2016)

Fig. 9

Plot of stand variables in each stand that contribute significantly to clustering

Fig. 10

The plot of studentized residuals against the predicted height for A. araucana. The studentized residuals were obtained by fitting model SM12 AA (a) without weighting, with weight ${w}_{i}=1/DB{H}_{i}^{k}$, where (b) k = 0.5, (c) k = 1, (d) k = 1.5, (e) k = 2, (f) k = 2.5, and (g) k = 3

Fig. 11

The plot of studentized residuals against the predicted height for A. araucana. The studentized residuals were obtained by fitting model Eq. 1 without weighting (left) and with weighting factor ${w}_{i}=1/DB{H}_{i}$ (right)

Fig. 12

The plot of studentized residuals against the predicted height for A. araucana. The studentized residuals were obtained by fitting the simple NLME H–D model for AA without weighting (left) and with power variance function (right)

Fig. 13

The plot of studentized residuals against the predicted height for A. araucana. The studentized residuals were obtained by fitting the generalized NLME H–D model for AA without weighting (left) and with power variance function (right)

Fig. 14

The plot of studentized residuals against the predicted height for N. pumilio. The studentized residuals were obtained by fitting model SM1 NP (a) without weighting, with weight ${w}_{i}=1/DB{H}_{i}^{k}$, where (b) k = 0.5, (c) k = 1, (d) k = 1.5, (e) k = 2, (f) k = 2.5, and (g) k = 3

Fig. 15

The plot of studentized residuals against the predicted height for N. pumilio. The studentized residuals were obtained by fitting model Eq. 3 without weighting (left) and with weighting factor ${w}_{i}=1/DB{H}_{i}$ (right)

Fig. 16

The plot of studentized residuals against the predicted height for N. pumilio. The studentized residuals were obtained by fitting the simple NLME H–D model for NP without weighting (left) and with power variance function (right)

Fig. 17

The plot of studentized residuals against the predicted height for N. pumilio. The studentized residuals were obtained by fitting the generalized NLME H–D model for NP without weighting (left) and with power variance function (right)