A percentile-based estimator for the log-logistic function: Application to forestry

The data used in this study were collected from two natural forest stands (multiple species composition with indeterminate age structure) and two monoculture stands of Gmelina arborea Roxb. and Tectona grandis Linn. f. One of the natural stands is situated in Ikrigon Forest Reserve, Cross River State, Nigeria and lies between Latitude 6°17.597′–6°17.862′ N and Longitude 8°35.597′–8°35.276′ E with an area of 542.7 ha. The second natural forest and part of the G. arborea stand are in Oluwa Forest Reserve, Ondo State, Nigeria. The natural forest lies between 6.83°–6.91° N and Longitude 4.52°–4.59° E with an area of 629 km². The G. arborea stand lies between Latitude 6°55′–7°20′ N and Longitude 3°45′–4°32′ E and occupies an area of 87,816 ha (Onyekwelu, 2001). T. grandis and G. arborea data were also collected from Omo Forest Reserve. Omo Forest Reserve is situated in Ogun State, Nigeria between Latitude 6°35′–7°05′ N and Longitude 4°10′–4°19′ E with an area of 130,500 ha.

Data were collected from 75 and 35 temporary sample plots (TSPs) in G. arborea and T. grandis stands, respectively; while 10 and 7 TSPs were used in Ikrigon and Oluwa natural forests, respectively. The plot sizes ranged from 0.04 to 0.25 ha, to achieve at least 30 trees in a plot. All living trees within the plots of the monoculture stands were measured, while in the natural forest, only trees with diameter at breast height (Dbh) ≥ 10 cm were measured. Details of the sampling procedures used for the data collection in Ikrigon natural forest, Oluwa natural forest, G. arborea and T. grandis can be found in T.E. Adeniyi (unpubl.) and Ogana et al. (2015, 2017), and Chukwu & Osho (2018), respectively. The descriptive statistics of the Dbh, total height (Ht), quadratic mean diameter (dg), dominant height (i.e., average height of the 100 tallest trees per ha, H_o), basal area per ha (G), density (N) and stand age (A) are presented in Table 1.

The probability density function (pdf) and cumulative distribution function (cdf) of the LogL are expressed as: (1)f(x)=αβ(xβ)α−1[1+(xβ)α]−2f\left(x \right) = {\alpha \over \beta}{\left({{x \over \beta}} \right)^{\alpha - 1}}{\left[ {1 + {{\left({{x \over \beta}} \right)}^\alpha}} \right]^{- 2}}(2)F(x)=[1+(βx)α]−1F\left(x \right) = {\left[ {1 + {{\left({{\beta \over x}} \right)}^\alpha}} \right]^{- 1}} where f(x) = pdf; F(x) = cdf α = shape parameters (α > 0); β = scale parameter (β > 0). The LogL has a closed-form cdf as to facilitate the estimation of tree proportion in different diameter classes without resorting to numerical integration.

The idea of the percentile estimator for the LogL stems from the principle developed for the Weibull distribution by Clutter et al. (1983). The idea is, if two sample percentiles are known, each can be related to its corresponding LogL cdf, and the resulting equations are iteratively solved to get the estimates of shape (a) and scale (β). Thus, the equations were developed as follows: “Let X_p represent the p-percentile value in the sample so that 100p-percentile of the sample values are less than X_p”. Therefore, from the definition of the LogL cdf: (3)P=[1+(βXp)α]−1P = {\left[ {1 + {{\left({{\beta \over {{X_p}}}} \right)}^\alpha}} \right]^{- 1}}Equation (3) was then solved for X_p to give equation (4)(4)Xp=β(1P−1)−1/α{X_p} = \beta {\left({{1 \over P} - 1} \right)^{- 1/\alpha}} The parameters of the LogL distribution were derived from equation (4). Since this is the first attempt of applying this percentile-based estimator for the LogL, several sample percentiles were considered. The following sample percentiles were analysed:

• -

  Method 1 (M1): 25th and 75th sample percentiles

• -

  Method 2 (M2): 17th and 97th sample percentiles

• -

  Method 3 (M3): 24th and 93rd sample percentiles

• -

  Method 4 (M4): 40th and 80th sample percentiles

• -

  Method 5 (M5): 30th and 70th sample percentiles

• -

  Method 6 (M6): 35th and 65th sample percentiles

• -

  Method 7 (M7): 40th and 60th sample percentiles

• -

  Method 8 (M8): 33rd and 67th sample percentiles

(5)X25=β(10.25−1)−1/α{X_{25}} = \beta {\left({{1 \over {0.25}} - 1} \right)^{- 1/\alpha}}

(6)X75=β(10.75−1)−1/α{X_{75}} = \beta {\left({{1 \over {0.75}} - 1} \right)^{- 1/\alpha}}

To ascertain the suitability of the percentile method, it was compared with the well-known maximum likelihood estimator (MLE). The MLE involves taking the partial derivative of the log-likelihood function of LogL with respect to the parameters (α and β) and setting the expression equal to zero. The resulting function is solved by numerical algorithm to get the estimates of the parameters. The log-likelihood function (ℓ) of the LogL is expressed as: (7)log(l)=nlog(α)−nlog(β)+(α−1)∑i=1nlog(xi)−n(α−1)log(β)−2∑i=1n[1+(xiβ)α]log \left(l \right) = nlog \left(\alpha \right) - nlog \left(\beta \right) + \left({\alpha - 1} \right)\sum\limits_{i = 1}^n {log \left({{x_i}} \right) - n\left({\alpha - 1} \right)log \left(\beta \right) - 2\sum\limits_{i = 1}^n {\left[ {1 + {{\left({{{{x_i}} \over \beta}} \right)}^\alpha}} \right]}} the partial derivative of equation (7) with respect to α and β will give equations (8) and (9), respectively: (8)∂l∂α=nα+∑i=1nlog(xi)−nlog(β)−2∑i=1n(xiβ)αlog(xiβ)1+(xiβ)α=0{{\partial l} \over {\partial \alpha}} = {n \over \alpha} + \sum\limits_{i = 1}^n {log \left({{x_i}} \right) - nlog \left(\beta \right) - 2\sum\limits_{i = 1}^n {{{{{\left({{{{x_i}} \over \beta}} \right)}^\alpha}log \left({{{{x_i}} \over \beta}} \right)} \over {1 + {{\left({{{{x_i}} \over \beta}} \right)}^\alpha}}} = 0}}(9)∂l∂β=−nβ−n(α−1)β+2αβ∑i=1n(xiβ)α1+(xiβ)α=0{{\partial l} \over {\partial \beta}} = - {n \over \beta} - {{n\left({\alpha - 1} \right)} \over \beta} + {{2\alpha} \over \beta}\sum\limits_{i = 1}^n {{{{{\left({{{{x_i}} \over \beta}} \right)}^\alpha}} \over {1 + {{\left({{{{x_i}} \over \beta}} \right)}^\alpha}}} = 0} where ∂ is the symbol of the partial derivative, ℓ represents the likelihood function, x_i is diameter at breast height and i ranged from 1 to n, other parameters are previously defined in equation (1). Since there are no explicit solutions to equations (8) and (9), the estimates were obtained numerically using the mledist function from the fitdistrplus package (Delignette-Muller & Dutang, 2015) implemented in R (R Core Team, 2017). The R code for the MLE is presented in Appendix 1.

Three goodness-of-fit statistics were used to evaluate the suitability of the percentile estimator. For each method, the Kolmogorov-Smirnov (KS), Anderson-Darling (AD) and Cramer-von Mises (W²) statistics were computed. The smaller the statistics are, the better the method.

KS statistic: (10)KS=max{max1≤i≤ni[Fn(xi)−F0(xj)],max1≤i≤ni[Fo(xj)−Fn(xi−1)]}KS = max \left\{{{{max}_{1 \le i \le {n_i}}}\left[ {{F_n}\left({{x_i}} \right) - {F_0}\left({{x_j}} \right)} \right],{{max}_{1 \le i \le {n_i}}}\left[ {{F_o}\left({{x_j}} \right) - {F_n}\left({{x_{i - 1}}} \right)} \right]} \right\}

AD statistic: (11)AD=−ni−∑j=1ni(2j−1)[ln(Fo(xj))+ln(1−Fn(xi−1))]/niAD = - {n_i} - \sum\limits_{j = 1}^{{n_i}} {\left({2j - 1} \right)\left[ {ln \left({{F_o}\left({{x_j}} \right)} \right) + ln \left({1 - {F_n}\left({{x_{i - 1}}} \right)} \right)} \right]/{n_i}}

W²statistic: (12)W2=∑i=1n{F^(xi)−i−0.5n}2+112n{W^2} = \sum\limits_{i = 1}^n {{{\left\{{\hat F\left({{x_i}} \right) - {{i - 0.5} \over n}} \right\}}^2} + {1 \over {12n}}} where F(x_i) = observed cumulative frequency distribution for x_i (i ranged from 1 to n), F₀(x_i) = theoretical cumulative frequency distribution.

The relative rank introduced by Poudel & Cao (2013) and recently used by Sun et al. (2019) was used for this study. It is expressed as: (13)Ri=1+(m−1)(Si−Smin)Smax−Smin{R_i} = 1 + {{\left({m - 1} \right)\left({{S_i} - {S_{min}}} \right)} \over {{S_{max}} - {S_{min}}}} where R_i = relative rank of method i (i = 1, 2, …, m); m = number of methods assessed (8 percentiles plus the MLE, i.e., 9 altogether), S_i = evaluation statistics value of method i; S_min and S_max = minimum and maximum value of S_i, respectively. The relative rank is a real number between 1 (best) and 9 (worst). For each method, the relative ranks were summed across the three goodness-of-fit statistics to ascertain the best percentile point for the estimator.

A scatterplot was used to establish the form of association between the dependent (LogL parameters and diameter percentiles) and independent variables (stand variables, e.g., quadratic mean diameter (dg), dominant height (H_o), stand density (N), basal area per ha (G), stand age (A), etc.). Thereafter, a stepwise linear regression technique was used in the modelling process. This ensured that only the best predictor variables were included in the models. Relating the parameters directly with stand variables is known as the parameter prediction method (PPM); while modelling the diameter percentiles with stand variables from which the LogL parameters could be derived is referred to as the parameter recovery model (PRM). Both PPM and PRM were evaluated to determine the best option for the LogL distribution. The models were of this form: (14)Yi=f(Xi,β)+εi{Y_i} = f\left({{X_i},\beta} \right) + {\varepsilon _i} where Y_i is a vector of the dependent variables (parameters of the LogL distribution and various diameter percentiles), X_i is the vector of the predictor variables (quadratic mean diameter, dominant height, stand density, basal area per ha, stand age, etc.) and is the error term.

Models were only developed for the G. arborea and T. grandis stands because a sufficient number of plots is required for effective modelling. There were more than 30 plots from these stands compared to the few plots in the natural stands (10 and 7 in Ikrigon and Oluwa, respectively). Since sufficient data are not available for modelling, a k-fold cross-validation was used to ascertain the predictive ability of the models (Sileshi, 2014). Different fit indices, such as root mean square error (RMSE), coefficient of variation (CV), corrected Akaike Information Criterion (AICc), Bayesian Information Criterion (BIC) and adjusted coefficient of determination ( R¯2{\bar R^2} ) were used to evaluate the cross-validation results.

The estimated parameters of the LogL from the different percentiles and those from the MLE are presented in Table 2 and Table 3 for the natural forests and plantations, respectively. Little variation was observed in the values of the parameters of the LogL estimated from the MLE and percentiles methods. Evaluation of the performance of the methods based on the fit statistics showed that MLE had the smallest relative rank sum (R) of 3.79 and 3.32 for Ikrigon and Oluwa natural forests, respectively. This was followed by method 4 (M4: 40th and 80th diameter percentiles), method 1 (M1: 25th and 75th) and method 7 (M7: 40th and 60th) with respect to Ikrigon forest. In Olwua natural forest, M4 had the largest R. The Cramer-von Mises (W²) statistics of M1, M4 and M7 were smaller relative to the value of the MLE in Ikrigon natural forest. Just as with KS and AD, W² statistic is an important index that gives the overall squared difference between empirical distribution and the theoretical cdf (Wang, 2005).

Similarly, in G. arborea stand, the MLE had the smallest relative rank sum of 3.05 which was closely followed by M4 with a value of 3.47. The performances of M1 and M7 were equally good. In the case of T. grandis stand, M4 had the smallest relative rank sum of 3.00 whereas MLE had 3.52. This was followed by M1 and M5 (30th and 70th diameter percentiles). Method 2 (17th and 97th diameter percentiles) with the largest relative rank sum had the worst performance in three out of the four stands. The 17th and 97th sample percentiles were recommended for the Weibull distribution by Shiver (1988). Thus, the study shows that these diameter percentiles are unsuitable for the LogL distribution. Gorgoso et al. (2007) also observed a poor fit with the percentiles for the Weibull distribution in the birch-dominated stands of Northwest Spain. Also, the most efficient sample percentiles – 24th and 93rd recommended for the Weibull function by Dubey (1967) performed poorly (M3) for the LogL function.

The graph of the relative frequency of trees against the equal diameter class of the best diameter percentiles and MLE for the four stands are shown in Figure 1a, b, c and d. The graph showed that the estimation methods produce a good representation of the four stands, that is, the reverse J-shaped and the Gaussian distribution typical of natural forests (uneven-aged) and plantations (even-aged), respectively. However, the relative frequency of trees was underestimated in the diameter class of < 20 cm in Oluwa natural forest (Figure 1b) and overestimated in the middle diameter class of G. arborea stand.

Figure 1

The results obtained with the LogL function in this study were better than those reported by Ogana & Wali (2018) for similar natural forest stands. Studies on LogL in forestry have used either the MLE (Ogana & Wali, 2018) or the least square regression method (Ogana & Dau, 2019). The least square method may seem handy, it is however, time-consuming. One important advantage of the MLE is the ability to evaluate the reliability of the estimates of the parameters using the standard error which can be derived from the hessian matrix. However, the MLE requires specifying the log-likelihood function for the distribution which may not converge due to poor initial guess for the parameters. The MLE is also affected by sample size (Shiver, 1988). Contrary to the MLE and least square method, the percentile estimator is more simplified, handier and does not require a complex iterative procedure. Therefore, when considering simplicity, vis-à-vis the ease of estimating the LogL parameters, the percentile estimator can be considered as a good alternative.

Two sets of models for predicting the diameter distributions of G. arborea and T. grandis stands were developed. The first set of models relate the best diameter percentiles with stand variables (parameter recovery model, PRM) in a simple linear relationship wherein the shape and scale parameters of the LogL were recovered. The quadratic mean diameter (dg) was the only suitable predictor variable for the dependent variables (diameter percentiles). The second set of models relate the LogL parameters directly with stand variables (parameter prediction model, PPM). The scale (β) parameter was predicted from the quadratic mean diameter (dg) in a simple linear relationship. The shape (α) parameter was predicted from the inverse of the natural logarithm of the ratio of the 25th diameter percentile and dg i.e., α=a+b1ln(P25/dg);\alpha = a + b{1 \over {ln \left({\raise0.7ex\hbox{${P25}$} \!\mathord{\left/ {\vphantom {{P25} {dg}}}\right.}\!\lower0.7ex\hbox{${dg}$}}\right)}}; where a and b are the estimated parameters. The parameters and standard errors of the models are presented in Table 4. The results of the 5-fold cross-validation showed that both model forms had a relatively low RMSE, CV, AICc, BIC and a high R¯2{\bar R^2} . Thus, they can be relied on for prediction purpose. Sun et al. (2019) also found dg as one of the most important predictors of the Weibull parameters in the mixed stand of oak and pine.

Table 4

Species	DV*	Parameter	Estimate	SE	5-fold Cross-validation
					RMSE	CV	AICc	BIC	R¯2{\bar R^2}
G. arborea	P25	b0	−0.590	1.570	4.830	4.801	116.937	123.460	0.533
		b1	0.632	0.069
	P40	b0	−2.975	1.676	5.539	5.441	126.491	133.015	0.677
		b1	0.912	0.074
	P60	b0	−0.768	1.159	2.519	2.615	72.663	79.186	0.854
		b1	1.050	0.051
	P75	b0	−0.677	1.116	2.446	2.418	67.149	73.672	0.897
		b1	1.232	0.049
	P80	b0	−0.468	1.148	2.567	2.561	71.250	77.774	0.900
		b1	1.288	0.051
	Shape	b0	1.952	0.124	0.138	0.147	−146.720	−140.196	0.764
		b1	−0.826	0.054
	Scale	b0	−0.753	0.955	1.803	1.768	44.331	50.855	0.860
		b1	0.888	0.042
T. grandis	P25	b0	−3.885	1.374	0.917	0.874	−1.105	2.556	0.844
		b1	0.975	0.074
	P40	b0	1.562	1.189	0.668	0.652	−10.687	−7.025	0.874
		b1	0.951	0.064
	P60	b0	0.263	1.016	0.495	0.475	−21.063	−17.401	0.912
		b1	0.995	0.055
	P75	b0	2.704	1.232	0.671	0.704	−8.334	−4.672	0.873
		b1	0.982	0.066
	P80	b0	4.946	1.331	0.786	0.829	−3.200	0.461	0.833
		b1	0.905	0.071
	Shape	b0	2.514	0.397	0.368	0.378	−30.438	−26.776	0.773
		b1	−1.028	0.098
	Scale	b0	−1.306	0.625	0.186	0.181	−53.142	−49.480	0.965
		b1	1.002	0.033

*

DV = dependent variables

The graphs of the application of the PRM and PPM to independent stands of G. arborea (0.08 ha) and T. grandis (0.125 ha) are shown in Figure 2a and b. The predicted diameter distributions were comparable to the observed diameter distributions of the stands. The first PRM (models from the 25th and 75th percentiles) had KS, AD and W² of 0.0670, 0.7583 and 0.0751, respectively in the 0.08 ha G. arborea stand. These values were respectively, 0.0778, 0.6879 and 0.0782 in the second PRM (models from the 40th and 80th percentiles). The PPM had 0.0727, 0.5968 and 0.0669, respectively. In the T. grandis, the KS, AD and W² of the PRM were 0.0993, 0.9301 and 0.1780, respectively. Whereas the PPM had 0.1178, 2.0522 and 0.3479, respectively. The second PRM could not predict the diameter distribution of the stand; as such, it was not represented in the graph.

Figure 2

The performances of the PRM and PPM were relatively the same for the LogL distribution. One major argument with the PPM is the difficulty of developing models that would account for high variations in the parameters of distributions (Borders & Patterson, 1990). While this statement may be valid for the shape parameters of some distributions, such as the well-known Weibull, it is not application to the LogL function. The models predicting the shape and scale parameters of the LogL explained more than 0.75 (> 75%) of the variation in the parameters for the stands. Thus, both the PRM and PPM could be used for the LogL in the G. arborea and T. grandis stands.