Introduction
Site quality is one of the most important factors influencing tree and forest mass growth and, therefore, stand productivity. According to the ^{Food and Agriculture Organization (FAO) (1994)}, this term is used to indicate the productivity of a site for a determined forest species (^{Clutter et al., 1983}). In effect, site classification is based on climate, edaphic and vegetation variables that condition growth and yield. Generally, indicators intrinsic to the particular stand are used, dominant height being the most appropriate due to its low dependency on density and habitual silviculture treatments and its wide correlation with total production in terms of timber volume (^{Savill et al., 1997}).
For the majority of research, the dominant heightage relationship is the most practical, consistent and useful measurement as a site index indicator. This finding coincides with ^{Fassola and Wabo (1993)}, who stated that dominant height is a frequently used parameter for indicating site index and which, in this particular case, was used through the development of a family of heightage curves, where the site index is the dominant height at a determined age.
^{Raulier et al. (2003)} suggested that site index curves could be derived from the heightage relationships when there are at least three sources of different data: temporal plots, stem analysis and permanent plots. Data obtained from permanent plots show exact pattern changes in the dominant height of a stand over time, taking into consideration tree replacement dynamics within the dominant strata (^{Raulier et al., 2003}). Dominant height intervals for the base age are divided into quality classes, taking into account that each quality class should have a range that is limited by height. In general, the range used is 2 or 3 m (^{Philip, 1994}). According to ^{Clutter et al. (1983)}, three methods have been used for determining site index: the guide curve, parametric prediction and the difference equation method.
Before the development of computerized methods, the most common technique to determine site index was the guide curve method (^{Savill et al., 1997}). This method incorporated elements of subjectivity that made it difficult to perform statistical tests on the curve's goodness of fit. These problems are easily overcome when using mathematical models. Currently, site index is examined by applying constantly evolving mathematical models.
The curves derived from the site index models should meet a series of properties, notably polymorphism, a sigmoid growth trend with an inflection point and the capacity to achieve a horizontal asymptotic at advanced ages; should present a logical response that is invariable regarding the simulation path and reference age; and should have a small number of parameters (^{Cieszewski and Bailey, 2000}). Models expressed in finite differences constitute the most used current method for fitting height growth equations, as they guarantee the vast majority of the requirements exacted from these functions. This method permits the estimation of the dominant height of a stand at a determined age based on the known dominant height of the stand at any other age. Selecting the parameter to be eliminated determines the behavior of the model: anamorphic or polymorphic curves. The finite differences method for developing the site index is described in more depth in ^{Clutter et al. (1983)} and ^{Torres and Magana (2001)}.
The majority of research regarding the use of finite differences for calculating the site index is focused on species of pine such as Pinus sylvestris L. (^{BravoOviedo et al., 2004}), Pinus pinaster Ait. (^{López and Valles, 2009}), Pinus durangensis M., Pinus pinea L. or Pinus uncinata Ram. (^{Calama et al., 2004}).
With regard to Tectona grandis, (teak), several studies on site index have been performed in countries as diverse as Tanzania (^{Malende and Temu, 1990}), northern Ghana (^{Nunifu and Murchison, 1999}), Costa Rica (^{Bermejo et al., 2004}) and India (^{Upadhyay et al., 2005}). In the state of Tabasco, Mexico, T. grandis is one of the most common plantation species. It is considered as a suitable species for the rapid production of large volumes of wood, firewood, posts and other products. However, at present, no equation has been developed to evaluate the productive capability of teak plantations in this Mexican state. The aim of this study was to determine the site index with finite difference equations using data from permanent plots of T. grandis in the state of Tabasco, Mexico.
Materials and methods
Study area
The study area consisted of 10 plantations of T. grandis within the municipalities of Cunduacán, Teapa, Jalapa, Cárdenas and Balancán in the state of Tabasco, southeastern Mexico (Figure 1). The names of the sampling locations were Rancho Hawail, Rancho la Reforma, Falcón 1000, Falcón 300, Bocanegra, Rancho San Agustín, Rancho Bellavista, C16 and Santandreu.
Tabasco is located between lat 18°39′ north and 17°15′ south and long 91°00′ east and 94°07′ west. The mean annual temperature is 26 °C, and the annual total rainfall can reach 4000 mm. The following forest types are present in Tabasco: high evergreen forest, medium evergreen forest, low deciduous forest, low flooded forest (canacohital), Logwood (Haematoxylum campechianum) and gallery forest. Regarding edaphology in the study area, several soil types are present, the most important being Vertisols, Regosols, Solonchak, Gleysols, Cambisols, Fluvisols, Rendzinas and Acrisols. Two of the most important rivers of Mexico run through Tabasco: the MezcalapaGrijalva and the Usumacinta; furthermore, Tabasco accounts for approximately 30% of the total surface water flow in Mexico.
Data
Data were obtained from 10 teak plantations with 35 permanent sampling plots. The oldest plantations measured were established in 1994 and the youngest in 1999. Measurements were taken in 2003, 2004, 2005 and 2006 during the months of March, April, June, September and October. The diameter at breast height (d) was recorded (cm) with calipers and the total height (h) with a Haga hypsometer (m). The majority of the measurements were taken in plots of 100 trees, from each individual tree within a 16 m^{2} framework, representing a total of 625 trees ha^{−1}. However, some plots varied in surface area, owing to different plantation frames, with the space between trees ranging from 5.75 to 36 m^{2} tree^{−1}.
The data were analyzed to eliminate any possible errors until achieving a depurated base of heightdiameter pairs. The main dasometric variables (dominant diameter, mean height, mean diameter) were calculated using the recorded dendrometric data. The total number of trees measured in all plantations was 7260 trees; the min. and max. values for diameter at breast height were 2 and 42 cm, with a standard deviation of 2.95; the min. and max. values of total tree height were 2 and 29 m, with a standard deviation of 1.31; and the dominant diameter values ranged from 10 to 34 cm, with a standard deviation of 2.24.
Statistical analysis
To develop the heightdiameter (hd) relationship 21 hd models were fitted (Table 1). The fit was implemented using the MODEL SAS/STAT® procedure with the Marquardt algorithm. The initial fitting parameters were those used by the respective authors in their corresponding studies.
Model  Equation  

Stage (1975) 

Eq[1] 
Meyer (1940) 

Eq.[2] 
Bates and Watts (1980) 

Eq.[3] 
Vanclay (1995) 

Eq.[4] 
Curtis (1967) 

Eq.[5] 
Curtis II (1967) 

Eq.[6] 
Prodan (1965) 

Eq.[7] 
Huang and Titus (1992) 

Eq.[8] 
Sibbesen (1981) 

Eq.[9] 
Curtis et al.(1981) 

Eq.[10] 
Tang (1994) 

Eq.[11] 
Huang and Titus II (1992) 

Eq.[12] 
Seber and Wild I (1989) 


Ratkowsky and Reedy (1986) 

Eq.[14] 
Weilbull and Bailey (1979) 

Eq.[15] 
Chapman and Richards (1959) 

Eq.[16] 
Zeide (1992) 

Eq.[17] 
Richards (1959 

Eq.[18] 
Bailey (1979) 

Eq.[19] 
Seber and Wild II (1989) 

Eq.[20] 
Seber and Wild III (1989) 

Eq.[21] 
Note. h: estimated height; d: normal observed data; βi: regression parameters.
Once the hd models were fitted for each plantation, those that presented the best R^{2} (coefficient of determination), R^{2}_{adj} (the adjusted coefficient of determination) and mean squared error (MSE) considered for each plantation were selected. Once the model expression was selected for each plantation, the Assmann estimated dominant height values were calculated from the dominant diameter data.
The method used to determine the site index curves involved the construction of multiple dominant heightage series based on the data obtained from the permanent forest plots. Thus, five equations (Table 2) were selected as function candidates for fitting the dominant height growth model. Four were based on the model of finite difference equations. The Schumacher function has been used to determine site index curves by means of the guide curve model (^{Clutter et al., 1983}) but can be transformed to obtain polymorphic curves. The Richards function has been widely used to determine forest growth (^{Amaro et al., 1998}). The Sloboda function is strictly increasing and is unusual in that it is nonlinear for parameters (^{Kiviste et al., 2002}). The ^{McDillAmateis (1992)} function does not take an integral form and directly expresses itself as a finite differences equation.
Original function  Free parameters  Differences equation 

(1) Schumacher(1930)  a 


b 

(2.1) Mc Dill Amateis (1992)  b 



(3.1) Richards(1959)  b 



(4.1) Sloboda(1971)  b 


Note: H1, H2 dominant heights at ages T1 and T2; a, b, c and d: function fitting parameters.
All of these possible functions have a common characteristic in that the predictions are all invariant over time, meaning that projections made over the same period of time are equivalent, without taking into account period length or the number of projection intervals (^{Palahí et al., 2004}).
A comparative analysis of R^{2}, MSE, and squared sum of errors (SSE) was performed to select the most suitable function found from the fitting of the five finite difference growth equations (1.1.1 to 4.1.1). The Proc NLIN procedure was used in the process (SAS/STAT 9.1. SAS Institute Inc., Cary, NC. USA.).
To validate the site index models, 30% of the observed values were randomly selected. Essentially, the aforementioned values for each of the finite difference equations were calculated using the global parameters that were obtained in the initial fitting of models. Subsequently, the following statistical parameters were determined: the absolute mean of the residuals (AMRES), the root mean squared error (RMSE), the coefficient of determination of the models (^{BravoOviedo et al., 2004}) and the model efficiency (EF) to select the best model. The expressions of these parameters are:
Awhere Y_{i}: observed value;
The graphs of the residual values were examined visually, as were the graphs of curves fitted for distinct site indexes, superimposed on the curves of the dominant heights that were observed over time. Visual inspection is essential for selecting the most precise model, as the curve profiles can differ considerably, even when the goodness of fit statistics are similar.
Finally, the points that the curves pass through were determined. Therefore, the extreme tree heights at a standard age were used to calculate the site index. The height intervals were defined and given a specific range. Furthermore, a comparative analysis of the site index curves with other curves proposed by ^{Upadhyay et al. (2005)} was conducted, employing the same methodology of superimposing data used in developing this study's model. The polymorphic equation used by this author was efficiently developed by ^{Elfving and Kiviste (1998)}, Eq. [22].
where
b = 345.217 and b_{2} = 0.77858
where H is the dominant height at base age (T_{0}) T_{1} initial age and T_{2} final age, d and r parameters to be calculated.
Results
The cloud of heightdiameter points used in the fitting of the hd models (Figure 2) presented good performance for the data pairs of the inventories of the plots installed. Table 3 presents the statistical values of fitting of the best hd equations in each plantation. In all models, the parameters were significant; however, in the most of the cases, the models showed a low R^{2}, less than 50%. The Stage model was the worst of the models in plantation 27, while the Huang and Titus I was the best with an R^{2}=0.6638, R^{2}adj=0.6631 and MSE=2.9459 m in plantation 23.
Plantation  Model  Estimator  Estimator value  R^{2}  R^{2}adj pond  MSE 

8  Huang and Titus II (1992)  b0 b1 b2 
24.33486**(1.5160) 2.861315**(0.3738) 0.086675**(0.0125) 
0.4208  0.4197  8.0124 
14  Prodan (1965)  b0 b1 
0.41311**(0.0452) 0.270655**(0.0150) 
0.4708  0.4696  9.7463 
16  Huang and Titus I (1992)  b0 b1 
0.977148**(0.0483) 0.225341**(0.0184) 
0.5292  0.5284  4.5654 
17  Meyer (1940)  b0 b1 
20.42288**(0.4607) 0.094425**(0.0055) 
0.5858  0.5849  5.6353 
23  Huang and Titus I (1992)  b0 b1 
1.089689**(0.0411) 0.175669**(0.00207) 
0.6638  0.6631  2.9459 
25  Stage (1975)  b0 b1 
2.139497**(0.1216) 0.696906**(0.0231) 
0.4984  0.4978  4.5593 
26  Vanclay (1995)  b0 b1 
0.032548**(0.00192) 0.593742**(0.0254) 
0.5976  0.5969  3.3701 
27  Stage (1975)  b0 b1 
4.603588**(0.2510) 0.385238**(0.0173) 
0.2612  0.2608  4.5468 
29  Stage (1975)  b0 b1 
2.011869**(0.1813) 0.600192**(0.0308) 
0.5864  0.5849  2.2147 
46  Meyer (1940)  b0 b1 
20.12308**(0.7243) 0.083068**(0.00709) 
0.3995  0.3987  8.5861 
Note: The data between parentheses represent the approximate standard error; level of significance of the model parameters
^{**}P<0.0001, R^{2} adj pond: fitted, weighted coefficient of determination; MSE: mean squared error.
The number of plots used in the fitting was 25, and 10 were used in the validation of the site index models for T. grandis. Thus, the mean, min., and max. values and the standard deviations of plot dominant height (16.95, 11.06, 22.43 m, 2.12) and plot age (8.86, 3.63, 12.39 yr, 1.97) of the fitting and validation (17.18, 11.60, 22.5, 2.22) and (8.24, 2.72, 12.39 yr, 2.02) are presented.
Fitting of the site index models
All of the parameters of the site index models were significant; however, the nonlinear function presented the best results of the statistical parameters (Table 4).
Difference model  Estimated parameter  SSE  MSE  R^{2}  Validation  

a  b  c  d  AMRES  RMSE  EF %  
Schumacher (1930)  –  4.0126 (17.4330)  0.0593 (0.2307)  –  17.8725  0.2628  0.9964  0.0155  0.0913  0.9236 
Schumacher F10(1930)  3.7185 (0.4857)  –  0.2985 (0.1597)  –  16.9822  0.2497  0.9984  0.1713  0.6419  0.9384 
McDillAmateis (1992)  33.2296 (9.4645)  –  0.5471 (0.1620)  –  17.1054  0.2515  0.9983  0.0148  0.0868  0.9374 
Richards 1.2 (1959)  24.8918 (2.8775)  –  0.3348 (0.0542)  –  17.3270  0.2548  0.9983  0.0146  0.6628  0.9343 
Sloboda F6 (1971)  3.5941 (0.4069)  –  1.0738 (1.6663)  0.2137 (0.2600)  16.7740  0.2504  0.9984  0.0137  0.0803  0.9377 
Note: The data between parentheses represent approximated standard error; R^{2}: coefficient of determinations. Significance level of the models P< 0.0001.
This function was derived from the Sloboda equation with a free b parameter, consistent with biological criteria. The Schumacher and Richards functions were those which presented the highest values for the SSE and the MSE and the lowest coefficient of determination (R^{2}), while the other three functions were similar regarding these particular aspects. The highest error was observed in the Schumacher function followed by the McDillAmateis function (Table 4).
Table 4 presents also the validation parameters of the models. Notably, all of the models exhibited good results. Considering the fit and the validation results, the Sloboda model was selected because the AMRES and RMSE values were closest to the ideal values (0).
This equation allows the prediction of dominant height as a function of plantation age. Taking into account the significance level of this model (P<0.0001), it was assumed that it is valid for the prediction of dominant height as a function of age. It is worth mentioning that in the case of the selected model; the mean squared error (MSE) was low, providing an adequate level of accuracy. The model's goodness of fit, expressed by the coefficient of determination (R^{2}), was a fundamental element in the selection model (R2=0.9984). The residual values versus the dominant heights and age (Figure 3) displayed accuracy and precision in the selected model. In Figure 3, it can be observed that the model was quite precise, indicating the inexistence of heteroscedasticity and the absence of outliers in the graph. The Sloboda differential function presented the best results in the fitting and validation processes. Therefore, the final expression of the site index model is shown in Eq. [23].
where H_{1} = dominant height at age T_{1}; H_{2=}dominant height at age T_{2}.
To determine the values the curves should pass through, the max. and min. heights at the base age of 10 yr were used, namely, 11.42 m min. and 25.22 m max.. Then, five site indexes were determined by dividing the data into five ranges and taking the central value of each range, varying from low, medium and high quality, considering that when the age (T_{2}) was equal to the base age (T_{1}), the dominant height (H_{2}) was equal to the site index (SI). The attained heights, depending on station quality, were biologically adequate. When the interval was divided into the five quality classes, the respective mean quality values 12, 15, 18, 21, and 24 were obtained.
The graphic representation of the Sloboda function (figure 4) forced to pass through the data pairs (10, 12), (10, 15), (10, 18), (10, 21), (10, 24) determined the curves for trees in the studied teak plantations in Tabasco. The growth function developed by ^{Upadhyay et al. (2005)} compared with our model (Sloboda) is shown in Eq. [24].
where
and H_{1} is the dominant height at the base age of 10 yr.
Discussion
All of the equations were fitted with data from permanent forest plots in the teak plantations at very early ages, between 2.72 and 12.39 yr. Nevertheless, this fit does not represent a serious problem, as teak is a fastgrowing species and in this particular case grows within a zone characterized by abundant precipitation and where edaphic properties are suitable for the normal development of this tree species.
The model used for the site index has been fitted using the finite difference equations method. Among the five evaluated difference equations of growth rate, the equation developed by Sloboda was selected, as it presented the best results for the different analyses conducted and met all of the required characteristics so that a model would be able to predict the heights of the analyzed plantations. The McDillAmateis function also presented suitable characteristics but was finally ruled out, as it presented statistical values (AMRES and RMSE) greater than the selected model and with a lower degree of efficiency. The equation derived from the Schumacher F10 function presented the highest efficiency (0.9384) of all of the adjusted models; however, this model was not selected, as it presented a negative AMRES value and a considerably greater RMSE than the selected model. Both the McDillAmateis and Schumacher F10 models were quite good in predicting the dominant height of T. grandis in Tabasco; however, overall, the equation derived from Sloboda presented superior results in fitting the equations. ^{Vargas et al. (2013)} suggests that the selected model has all of the desirable characteristics when estimating station quality, given that it has already been applied with the same aim in other timber species.
The method used for the construction of the site index curves was selected by the function of the representativeness of the observed value. The estimated polymorphic curves of the site index reflect the tendency of the cloud dispersion interval of the observed values, suggesting that the proposed method is good for achieving the established aim. The site quality curves obtained in this study clearly demonstrate that Tabasco is an excellent site for the establishment of teak plantations. In addition, it is important to analyze the residual values, with the aim of studying diverse relationships within and outside the range of data. Therefore, the residuals should be represented graphically as functions of age and dominant height to visualize these distributions over time. The residual analysis undertaken in this study led to the assumption that in all cases, the selected model is fairly good as an estimator for predicting the site index.
In other studies performed in teak plantations in the department of Cordoba, Colombia, ^{Torres et al. (2012)} found similar results using the polymorphic variant of the Korf function, this model being superior to the others included in their study. The results presented high variability regarding site quality. In the specific case of the plantations of Cordoba, located in the Caribbean region of Colombia, the authors used 12 yr as the typical age, as they were young plantations consisting of trees between 3 and 22 yr old. Five site qualities attained a mean dominant height of 24.8 at the best sites, between 15.8 to 18.8 m at the intermediate sites and 9.8 m at the worst sites. The results found in Cordoba were similar to those found in the present study.
^{Mora and Meza (2004)} selected the generalized Schumacher model as the best for representing site index curves for teak plantations on the pacific coast of Costa Rica from a series of lineal and nonlineal fixed effect models, using data from permanent plots and stem analysis. The selected guide curve suggests initial rapid upward growth that remains relatively high without reaching growth stabilization, even at 40 yr of age. ^{JerezRico et al. (2011)} obtained a guide curve in Venezuela that coincides with the rapid initial growth of teak. The author generated a family of anamorphic curves based on mixed models from the Schumacher model in its lineal and nonlineal forms for site indexes 27, 24, 21, 18 and 15 m at the base age of 16 yr. The results analysis demonstrated a better fit for the nonlineal mixed models when compared with the other models in terms of bias and precision, demonstrating the convenience of using models that consider repeated measures in the same plot. This finding would indicate that no differences exist between the sites in Venezuela, Costa Rica and Tabasco, Mexico.
However, ^{Sajjaduzzaman et al. (2005)} fitted the Chapman Richard growth model for teak in Bangladesh, using 40yras the standard age and max. and min. dominant heights of 23 and 3.5 m, respectively. Thus, the results of the present study were superior to those found in Bangladesh. ^{Dupuy et al. (1999)} developed curves for the SudanGuianese savannah in Ivory Coast, generated from a sample of 200 permanent and temporal plots, with a max. age of 61 yr. The curves presented good growth, particularly during the initial yr. However, at approximately 10 yr, there was a slight reduction in the slopes of the curves belonging to the three smallest sites in Ivory Coast, indicating that from this age onward, these curves reached a level inferior to the equivalent curves in Tabasco, Mexico.
A comparison has been made between the Sloboda equation fitted in this study and the model developed by ^{Upadhyay et al. (2005)}. The Sloboda model presents a better data fit at ages less than the base age when compared with the curves produced by ^{Upadhyay et al. (2005)} using equation [24]. However, regarding ages greater than the base age, the curves of equation [24] tend to overestimate the dominant height with respect to the Sloboda model, resulting in the quality curves of equation [24] attaining the asymptote of the dominant height at a greater age when compared with the model developed in the present study. The overestimation of the dominant height in the ^{Upadhyay et al. (2005)} model, using the base age, could result in erroneous classifications when estimating the evolution of dominant height with age and, consequently, in site productivity calculations when applying this model to teak plantations in Tabasco.
A model is often considered suitable when the goodness of fit statistics coincide, are without bias, and present high efficiency throughout their processing. However, a good model should also be validated using data not used in the fitting. Even though the group of functions evaluated in this study present low bias, low error and high efficiency, 30% of the data from the studied plots are reserved for this purpose. The residual values were small, resulting in acceptable degrees of error for both short and longterm growth projections of teak trees in the state of Tabasco, Mexico. However, the site index study should occur at the end of the first rotation, with the aim of verifying whether the Sloboda model is still the most appropriate for these plantations.
In summary, the Sloboda model, considering the parameter b as free, has been demonstrated to be the most appropriate in both the process of fitting and results validation. This model illustrates the evolution with age of the dominant height of teak trees in the state of Tabasco, Mexico, taking into account aspects related to the growth habits of this particular species. Through this model, five site index curves have been defined (all with dominant heights between 12 and 24 m at a reference age of 10 yr), which will provide great assistance in classifying the productivity of new teak plantations in Tabasco. Therefore, with more reliable information regarding the dynamics of this species at the established sites, longterm monitoring will be indispensable, given that in the majority of cases, the high demand for timber products prevents trees reaching the established age for plantation rotation as a fundamental part of the integral management of these plantations.