Researchers realize that volume, the primary varia­ble that forest managers are interested in, is heavily de­pendent on both tree diameter and height. A traditional practice is to fit a marginal distribution to the diameter frequency data and then use an empirical height-diameter relationship to estimate the average height per diameter class and hence the volume ( Clutter & Al­li­son, 1974 ). Although this approach is satis­fac­to­­­ry in many situations, it is often not appropriate because the approach tends to ignore the natural relationships between tree diameters and heights by treating them separately ( Schreuder & Hafley, 1977 ; Tewari & von Gadow, 1999 ).

The traditional method does not quantify the distribution of heights for a given diameter and one approach for modelling the conditional height distribution for the different diameters is to use the height residuals ( Gaffrey, 1996 ). However, it is very seldom that the height residuals are homoscedastic and normally distributed. In most forest stands the variance about the diameter-height regression is not homogeneous ( Zucchini et al ., 2001 ).

One of the most important elements of forest structu­re is the relationship between tree diameters and heigh­ts because information about size-class distributions of the trees within a forest stand is important for estimating product yields ( Zucchini et al ., 2001 ). The size-class distribution influences the growth potential and hence the current and future economic value of a forest stand ( Knoebel & Burkhart, 1991 ). On the other hand, the social status of a tree, which reflects its further development in the aspects of growth and mor­ta­­li­ty, depends not only on its relative diameter, but also on its relative height in a stand ( Siipilehto, 2000 ). Furthermore, knowledge of the height variation, both between and within diameter classes, improves the chances of successfully imitating different types of thinnings ( Hafley & Buford, 1985 ).

For many years, the bivariate extension of the S _B dis­tribution, the S _BB ( Johnson, 1949a and 1949b ), has been the most commonly bivariate distribution used for modeling bivariate tree diameter-height frequency data (e.g. Schreuder & Hafley, 1977 ; Hafley & Buford, 1985 ; Knoebel & Burkhart, 1991 ; Siipileh­to, 1996 ; Tewari & von Gadow, 1997 ; Tewari & von Gadow, 1999 ; Zucchini et al ., 2001 ; Li et al ., 2002 ; Wang & Rennolls, 2007 , Gorgoso-Varela et al ., 2016 ; Ogana et al ., 2018 ; Ogana, 2018a ). The Johnson's S _BB is developed by applying a four-parameter logistic transformation to each of the component variables of a standard bivariate normal distribution ( Johnson, 1949b ; Wang, 2005 ; Gorgoso-Varela et al ., 2016 ). The fitting of the joint distribution of diameters and heights is called bivariate distribution modeling ( Ogana et al ., 2018 ). Bivariate distribution modelling allows for the generation of bivariate frequencies of tree diameter and height ( Tewari & von Gadow, 1997 ).

The parameters of the Johnson's S _BB distribution ha­ve been estimated with both the indirect and direct me­thods. The indirect method fits the bivariate Johnson's S _BB distribution by two-stage where the marginals are first fitted separately for the diameter and height; the estimates are then used to compute the parameter of association. Examples of the indirect method that have been reported in forestry literature include Conditional Maximum Likelihood (CML), moments, mode and Knoebel and Burkhart methods ( Hafley & Buford, 1985 ; Tewari & von Gadow, 1997 ; Knoebel & Burkhart, 1991; Siipilehto, 1996 ; Tewari & von Gadow, 1999 ; Cas­tedo-Dorado et al ., 2001; Li et al ., 2002 ; Oga­na, 2018a ). Of these indirect methods, CML and moments ha­ve been reported as the best methods ( Ogana, 2018a ). In the direct method, all the parameters in the bivaria­te dis­tri­bu­tion are estimated simultaneously ( Wang, 2005 ; Mønness, 2015 ; Gorgoso-Varela et al ., 2016 ; Oga­na et al ., 2018 ). However, these two methods (i.e., direct and indirect) of estimating the parameters of Johnson's S _BB distribution have never been compared using the same sample. The method used to calibrate the distribution matters in order to evaluate the accuracy of diameter and height predictions ( Ogana, 2018a ). Thus, the objective of the study was to com­pa­re the indirect fitting method for the Johnson's S _B based on two-stage method by the Conditional Maximum Likelihood and moments methods with the direct method where the parameters are estimated simultaneously.

The descriptive statistics i.e., the mean, maximum, minimum and standard deviation of the estimated parameters of the bivariate Johnson's S _BB distribution are presented in Table 2. The location parameter of the marginal height distribution was set to 1.3 for both direct and indirect methods across the three species. The shape, scale and rho (ρ) parameters of the indirect methods (CML and moments) show little variations. In fact, the standard deviations of the parameters' values were the same at least up to 1 decimal place. This indicates that both CML and moments are closely related. However, estimates from the direct method were different. The standard deviations of the estimated parameters are smaller than that of CML and moments of the indirect method.

The assessment of the suitability of indirect (CML and moments) and direct methods for tree height pre­diction showed that the indirect method, especially the CML method, had the highest R ² and the lowest RMSE, MEF, BIC and HQC in E. globulus, P. pinas­ter and P. radiata (shown in Table 3). The R ² , RMSE, MEF, BIC and HQC of CML ranged from 0.994 - 0.820, 1.454 - 1.676, 0.055 - 0.179, 4559 - 13400 and 4492 - 13335, respectively. However, the direct method had relatively lower R ² and larger RMSE, MEF, BIC and HQC in E. globulus, P. pinaster and P. radiata relative to the indirect methods. The number of plot in which the observed and predicted values were significant at 5% level (t-test) was 1 out of 308 plots (equivalent to 0.32%) in E. globulus for CML. Eight (2.59%) and 14 (4.45%) plots were significant for moments and the direct methods, respectively in E. globulus . There was no significant difference between the observed and predicted tree height by CML and moments across the 184 and 96 plots in P. pinaster and P. radiata , respectively. But 4 and 5 plots were significant with the direct method in P. pinaster and P. radiata , respecti­­ve­­ly at the specified level.

In the case of tree volume prediction, indirect (CML and moments) and direct methods performed relative­ly the same (Table 4). The methods had the same R ² , RMSE and MEF values up to 2 decimal places. The direct method had lowest BIC and HQC in E. globulus and P. radiata . The number of plot in which the obser­ved and predicted values was significant at 5% level for CML was like the result in tree height prediction (i.e., 1, 0, 0). However, 6, 6, and 3 plots were significant in E. globulus, P. pinaster and P. radiata , respectively for the direct method.

The mean of the residuals from the height-diame­ter median regression was computed across DBH cla­sses and plotted accordingly for the three species. The data were group into DBH classes of 5 cm interval and the mean residual prediction for each class by species was assessed. The graphs showed that CML and moments over- and under-estimated the tree height in the larger diameter classes (> 42.5 cm) in E. globulus and P. pinaster (Fig. 1a to c). Nevertheless, their prediction errors lie within the range of - 0.5 to 0.5, except in the larger DBH classes. Conversely, the direct method did not only over- and under-estimate tree height in the larger diameter classes but also over-estimate tree heights in the lower diameter classes in E. globulus . The mean residual height prediction for both direct and indirect method were relatively poor in P. radiata .

Furthermore, the graph of residual against predict­ed tree volume showed that CML and moments and the direct methods occupied the same horizontal band (Fig. 2). Their mean residuals plot lied within the range of -0.6 to 0.6, -0.6 to 0.6 and -0.4 to 0.2 in E. globulus, P. pinaster and P. randiata , respectively. Such plot is used to assess whether the residual is homoscedastic (i.e. constant variance) or heteroscedastic. The residual of the methods assessed were relatively constant.

The direct and indirect methods of fitting Johnson's S _BB have been evaluated. The methods follow similar trend across E. globulus, P. pinaster, and P. radiata with CML being the best for tree height and volume predictions. The number of plots rejected by t-test increased from CML to moments and to the direct method. The percentage of rejections in E. globulus was higher than the other species because some plots of coppice stands were considered, in which the minimum diameter considered was less than 5 cm (value of the minimum diameter considered for Pinus pinaster and Pinus radiata stands). There were also biases in the larger diameter classes (> 42.5 cm) for all methods. Kalbi et al . (2017) attributed this bias to the small number of trees in the larger diameter class. Few trees were found in the diameter classes > 42 cm across the three species. The Johnson's height-diameter function often perform well if parameter Φ i.e., the determinant of the regression curve is greater than one ( Tewari & von Gadow, 1999 ). In this study, Φ had values greater than one in some of the sample plots for CML, moments and direct methods. The model performed relatively well in these plots.

Furthermore, the values of the BIC and HQC diffe­rence show that the method of CML fits the diameter-height data better than direct methods for height prediction in E. globulus, P. pinaster and P. radiata. As a rule of thumb, a minimum ∆AIC/∆BIC/∆HQC of ≤ 2 is required for two models to be indistinguishable ( Ogana, 2018b ; Tewari & Singh, 2018). Though, the measures of the accuracy of the estimates (RMSE) and the relative measures of the model performance (MEF) look quite similar for all methods, the CML seems to be slightly better. Parallel observation was reported in Ogana (2018a) , who evaluated the performance of CML, moments, Knoebel and Burkhart (KB) and mode methods for height and volume predictions. The author found CML and moment to be the best methods for predicting tree height and volume in E. camaldulensis Dehn.

It is obvious from this study that when a complex fitting approach does not outperform another simpler one it is convenient to take the latter. Estimating the parameters of Johnson's S _BB by CML or moments is much easier and pose little computational difficulty. The CML and moments estimation procedure can even be carried out on Microsoft excel platform. However, the direct method involves complicated algorithm using maximum likelihood technique. This often requires writing the log-likelihood function and the specification of initial values for the parameters which may not achieve convergence ( Wang, 2005 ). Given the complexity of the likelihood function, its gradient is not tractable and then numerical optimization methods have to be used (i.e. through the 'optim' function), which may be less efficient than those based on derivatives. The computation procedure gets more complex as the number of parameters in the distribution increases. Bivariate Johnson's S _BB has nine parameters which makes fitting somewhat complicated. Often constraints are imposed on the two location and scale parameters of the bivariate distribution to make the estimates more plausible ( Wang & Rennolls, 2007 ). Mønness (2015) used maximum likelihood estimation to fit the Johnson's S _BB distribution and compared the result with power-normal and hyperbolic height prediction model. The prediction from power-normal and hyperbolic models were better than bivariate Johnson S _BB fitted with maximum likelihood. Other studies that have used maximum likelihood technique to fit the bivaria­te Johnson's S _BB distributions include Wang et al . (2008) , Gorgoso-Varela et al . (2016) and Ogana et al . (2018) . Satisfactory results were obtained by the authors for the Johnson's S _BB relative to other bivariate functions evaluated in their study.

The bivariate Johnson's S _BB height-diameter median regression remains the most applied function amidst several bivariate distributions in forestry literature. This is because of its flexibility and, more importantly, because its implied relationship between height and diameter is biologically reasonable ( Burkhart & Tome, 2012 ). The closed-form median regression of the Johnson's S _BB height-diameter model has also enhan­ced its continuous application to forestry. This median regression is frequently used for computing percentile lines wherein bounds on height are set. These percentile lines can be used to show how the variation in height decreases with increasing tree size for a specified diameter ( Tewari & von Gadow, 1999 ). The scatter plots of the Johnson's S _BB me­dian regression using indirect method fitted by CML from some sample plots in E. globulus, P. pinaster and P. radiata are presen­ted in Figures S1 to S3 [suppl.].

Modelling the joint distribution of diameter and height by Johnson S _BB remains an important tool for assessing the variation of tree height for a given diameter, detailed stand structure and volume estimation. Its accuracy is affected by the method of estimating the parameters of the distribution. In this study, we found the indirect (two-stage) method, especially by conditional maximum likelihood, to be the most suitable method for the bivariate Johnson's S _BB distribution. This method predicted tree height and volume in E. globulus, P. pinaster and P. radiata better than the direct method.