Sileshi F. Melesse · Temesgen Zewotir

Abstract The effect of tree age and climatic variables on stem radial growth of two hybrid clones of Eucalyptus was determined using longitudinal data from eastern South Africa. The stem radius of was measured weekly as the response variable. In addition to tree age, average weekly temperature, solar radiation, relative humidity and wind speed were simultaneously recorded with total rainfall at the site. An additive mixed effects model that incorporates a non-parametric smooth function was used. The results of the analysis indicate that the relationship between stem radius and each of the covariates can be explained by nonlinear functions. Models that account for the effect of clone and season together with their interaction in the parametric part of the additive mixed model were also fitted. The interaction between clone and season was not significant in all cases. For analyzing the joint effect all the covariates, additive mixed models that included two or more covariates were fitted. A significant effect of tree age was found in all cases. Although tree age was the key determinant of stem radial growth, weather variables also had a significant effect that was dependent on season.

Keywords Additive mixed effects · Dendrometer trial ·Parametric modelling · Penalized splines · Weather variables

Introduction

Eucalyptus trees, notable, economically valuable members of the Myrtaceae family (Kulheim et al. 2011), are among the main global sources of forest products such as fuelwood, transmission poles, plywood timber, pulp, and building materials (Melesse and Zewotir 2017a). Essential oils extracted from their leaves also have powerful medicinal properties. As a result, Eucalyptus comprises more than 500 species and is the most-planted genus of trees in the world (Demel 2000). Their fast growth makes them especially profitable for plantations to meet global demand for wood products, pulp and paper and enables harvesting within a short time. They can also resprout after harvesting for additional harvests. Because of this versatile nature of Eucalyptus species, many countries promote planting them (Melesse and Zewotir 2017a).

By the end of 20th century, Eucalyptus had become the most extensively planted hardwood species in the world(Turnbull 1999; Melesse and Zewotir 2017a). Eucalyptus is the most commonly planted genus in Africa, covering 22.4%of all planted area (FAO 2001; Melesse and Zewotir 2017a).Eucalyptus plantations in Africa are owned by small farmers for subsistence living and rural livelihoods in terms of fuel,poles and building materials and by large conglomerates for industrial wood. Industrial Eucalyptus plantations are 52%publicly owned, 34% privately owned and 14% other or unspecified, while of the non-industrial plantations, 62% are publicly owned, 9% privately owned and 29% other or not specified (FAO 2001; Melesse and Zewotir 2017a).Although 72% of the forest plantations are owned privately in South Africa (CIFOR 2000; Melesse and Zewotir 2017a),privately owned plantations are not well developed in the rest of Africa. However, some are being developed in countries such as Uganda, Mozambique, Tanzania and Rwanda. Eucalyptus trees are economically very important for pulp production and solid wood applications and a chief source of income for South Africa.

Eucalyptus grandis is the most widely commercially grown species in South Africa, making up about 73.8% of the total commercial forestry (Owen 2000). It is fast growing and has well-studied wood properties, but it is less resistance to drought and cold (Stanger et al. 2011) than Eucalyptus camaldulensis, which grows slowly. The rapid demand for hardwood, coupled with scarcity of land to grow E. grandis in South Africa has necessitated the expansion of plantations into colder sites where E. grandis might not survive. To solve this problem, breeders have crossed E. grandis with Eucalyptus camaldulensis and Eucalyptus urophylla to produce hybrids that are more adaptable and perform better than E. grandis in poor sites (Malan 1995; Melesse and Zewotir 2017b; Van Wyk et al. 1988).

Modelling juvenile tree growth is important in forest management to determine timber yield and long-term response of forest structure and dynamics to selective logging (Chazdon et al. 2010; Finegan et al. 1999; Herault et al.2010). The greatest potential for improving growth rates is during juvenile development (Watt et al. 2004), so appropriate juvenile development modelling is crucial for simulation models (Gang et al. 2011). Among studies that have considered modelling the growth of juvenile Eucalyptus tree are those of Melesse and Zewotir(2013a, b, 2015, 2017a, b, c). In studies (Farrow et al. 1994;Honkanen et al. 1999; Eyles et al. 2009; Turnbull et al. 2007)that used small plants or juvenile trees (aged ＜3 years),stem radius data was fitted using parametric regression methods for longitudinal data. However, none of these studies used nonparametric or semi-parametric approaches.Although parametric models provide a powerful tool for modelling the relationship between the responses and the covariates, they suffer from inflexibility when used to model complicated relationships between the response variables and covariates in various longitudinal studies (Lin and Carroll 2008). In parametric methods, the form of the underlying relationship must be known in advance except for the values of a finite number of parameters. That means the relationship between the mean of the longitudinal response and the covariates is fully parametric. In such situations, more flexible models that can adequately capture the pattern characterized by a function that is nonlinear in its parameters must be considered (Fitzmaurice et al. 2008).

The main drawback of parametric modelling is that it may be too restrictive or limited for many practical cases. This limitation has motivated a demand for developing nonparametric regression methods for analysis of longitudinal data.These methods can help to estimate a more flexible functional form between the responses and the covariates in the data.Consequently, complicated relationships between longitudinal responses and covariates can possibly be captured from the data. The main idea behind the nonparametric approach is to let the data ‘‘decide’ the most suitable form of the functions.Nonparametric and parametric regression methods should not be regarded as competitors, instead they complement each other (Wu and Zhang 2006). In some situations, nonparametric techniques can be used to validate or suggest a parametric model. A combination of both nonparametric and parametric methods is more powerful than any single method in many practical applications.

Although parametric models may be restrictive for some applications, nonparametric models may be too flexible to make concise conclusions in comparison with parsimonious parametric models. Semi-parametric models are good compromises and retain valuable features of both the parametric and nonparametric models (Fan and Li 2004).Some of these studies (Farrow et al. 1994; Melesse and Zewotir 2013a, b, 2015, 2017a, b, c) considered the longitudinal nature of the data. However, they all used parametric modelling; none used nonparametric or semiparametric approaches. The current study thus used a semiparametric approach to analyze the effect of climatic factors on the stem radial growth of hybrid clones of E.grandis × E. urophylla and E. grandis × E. camaldulenis.

Materials and methods

Data

A dendrometer trial focused on the growth of Eucalyptus grandis × E. urophylla (GU) and an E. gradis × E.camaldulensis (GC) hybrid clones established on Sappi landholdings near the town of KwaMbonambi in the KwaZulu-Natal Province of South Africa (Melesse and Zewotir 2013a, b, 2015, 2017a, b, c). The trial was designed to run over several years with separate growth monitoring phases. Nine trees were selected from each clone for intensive monitoring of radial growth (Drew 2004; Drew et al. 2009; Melesse and Zewotir 2013a, b, 2015, 2017a, b, c). Several wood characteristics,including stem radius measured with a caliper, for 18 trees,nine from each clone, were collected from April 2002 when the trees were 39 weeks old until August 2003 when trees were 107 weeks old. Data were compiled in a database for the entire trial.

From the 18 sampled trees (nine per clone), longitudinal data for 1242 weekly stem radial measurements were obtained. The response variable investigated in this study was the weekly stem radial growth, which is of interest because it can be used to understand the underlying processes of fiber development in fast-growing Eucalyptus plantations. In addition, the study of young trees may be very important in the selection of more-productive tree species. Refer to studies by Drew et al. (2009) and Melesse and Zewotir (2013a, b, 2015, 2017a, b, c) for details on data collection, field preparation and soil surveys.

Methods

The classical linear regression model and cross-sectional study only deal with the average change and does not provide information on how individuals change over time.Changes at the group and individual levels can be discovered only within the framework of a longitudinal study.

A unique feature of longitudinal data is that observations within the same individual are correlated. Applying methods that account for the effect of these correlations is vital to avoid an erroneous estimation of the variability of parameter estimates. The interdependence among measurements of the same individual can be captured in the mixed modelling framework.

The current data set consisted of repeated measurements of the same subjects over time. Numerous linear and nonlinear mixed effects have been proposed in the analysis of the longitudinal data (Verbeke and Molenberghs 1997, 2000; Fitzmaurice et al. 2004; Meng and Huang 2010). Models for the analysis of such data recognize the relationship between serial observations on the same unit.Most work on methods of repeated measures data has focused on data that can be modeled by an expectation function that is either linear or non-linear in its parameters(Laird and Ware 1982; Melesse and Zewotir 2015).

The additive model can be formulated by introducing the smoothing function of some predictor variables in the classical linear regression model.

where y is n × 1 vector of the response variable, X*is n × (p + 1) model matrix for the parametric components of the model, α is the corresponding (p + 1)× 1 parameter vector and the fj(·) is a smooth arbitrary function of a covariate xjand has the same dimension as y, ε is the n × 1 vector of random errors. The assumptions of the additive model are the same as the assumptions in the linear model except for the assumption of linearity. The error variance is constant, the error is normally distributed,and the errors are uncorrelated.

Additive mixed models

The insertion of the random effects into the additive model(1) gives us the additive mixed model.

where Z is the n × (q + 1) design matrix for random effects b with dimension (q + 1), ε is a vector of random error that is independent of b, ε ～N(0， R) and b ～N(0， Gθ). Both covariance matrices R and Gθare positive definite. These matrices are also assumed to depend on a parsimonious set of covariance parameters.The additive mixed model (AMM) that can have a nonnormal response is a generalized additive mixed model(GAMM). The GAMM has the following structure:

where G (·) is a monotonic differentiable function. The GAMM represents a model with higher flexibility and complexity, where mixed effects, smooth terms and nonnormal responses are included (Lin and Zhang 1999).These models can be viewed as additive extensions of the generalized linear mixed models.

Inference in generalized additive mixed models

Statistical inference in generalized additive mixed models includes estimations of the non-parametric functions fj(·), the smoothing parameters, λ, and all the variance components. In the case of Gaussian response and identity link function, the estimation of nonparametric functions, smoothing and variance parameters in the context of GAMM is achieved using restricted maximum likelihood (REML).

For a non-Gaussian response, penalized quasi likelihood(PQL; Breslow and Clayton 1993) and double penalized quasi likelihood (DPQL) are used to estimate the parameters and nonparametric function (Lin and Zhang 1999).Both PQL and DPQL originate from the maximum likelihood (ML) technique. The ML has direct application only in fixed models with a Gaussian response. The ML approach is also used in linear mixed models; however, the maximum likelihood estimators (MLEs) of variance are biased in general. As a result, restricted maximum likelihood (REML) estimators are used instead of MLEs.

Both ML and REML assume that the response is normally distributed. The assumption of normality is often easily violated in practice, making the likelihood inference difficult. In the absence of the random effects and errors distributions, the likelihood function is not available. Even in the presence of non-normal distributions of the random effects and errors with some unknown parameters, the likelihood function can be quite difficult to calculate and may not have appear analytical. Moreover, the distributional assumptions for any non-normal distribution may not hold in practice. These problems have led to the development of methods other than maximum likelihood, for example, quasi-likelihood, also known as Gaussian likelihood approach, which avoids the computational difficulty of the maximum likelihood method. Because the REML estimates can be derived from the quasi-likelihood (Heyde 1994), the Gaussian REML estimation can be considered as a quasi-likelihood method.

Penalized quasi-likelihood

When the exact likelihood function is computationally intractable, there are no simple solutions to get the parameter estimates. One possible option is to use numerical integration techniques. Some of these are Gaussian quadrature, numerical integration like Markov chain, Monte Carlo algorithms, stochastic approximations algorithms and penalized quasi-likelihood (PQL; Zuur et al. 2009). Penalized likelihood estimation has been proposed as a computationally simple alternative to methods based on numerical quadrature, especially when the number of random effects is relatively large (Fitzmaurice et al. 2004). By employing Laplace approximation, an approximated maximum likelihood, called PQL,can be obtained instead of exact likelihood. Penalized likelihood is essential in the case of non-normal data. (For details of the key concept in quasi-likelihood and Laplace’s approximation, refer to Lin and Zhang 1999; Valeria 2011;Zuur et al. 2009; Fitzmaurice et al. 2004 and Green and Silverman 1994).

Although several R packages (R Core Team 2017) are developed to fit GAMM, the most versatile that can handle modelling the correlation structure is the package mgcv(Wood 2006), which uses the nlme implementation of nonlinear mixed models. It also fits non-Gaussian responses by calling MASS’s generalized linear mixed model penalized quasi-likelihood (glmmPQL). The main advantage of this package is that it is possible to include serial and/or spatial correlation structures of the random effects.In this study, the package mgcv (Wood 2006) was used to fit the additive mixed models.

Results and discussion

Stem radius was measured every week for 69 weeks for 18 trees (nine from each clone) for a total of 1242 weekly measurements (9 × 69 = 621 measurements for GU clones; 621 measurements for GC clone were obtained).

The summary measures for the response variable, stem radius, and the climatic variables are presented in Table 1 and Table 2, respectively. From the summary measures,the mean stem radius for GU clone was larger than the mean stem radius for GC clone. The maximum stem radius for GU was higher than the maximum stem radius for GC clone. At the study area, the minimum temperature was 11.53 °C, and the maximum was 28.74 °C. Total rainfall ranged from 0 to 72 mm (Table 2).

Table 3 shows the results of the fitted additive mixed effects model (AMM) to each of the covariates used in the study. Both additive mixed model output with smoothed nonparametric part and linear model without smoothed part are presented for comparison. The nonlinearity was tested using a formal test by comparing a model specifying the smooth term with a model specifying a linear trend only.The difference between the two models (linear trend versus smooth terms) is tested using the likelihood ratio test statistics. As shown in Table 3, for all comparisons the pvalues were (p ＜0.005) very small. That means the relationship between stem radius and each covariate is nonlinear.

Moreover, by looking at the last column of Table 3, we can see the effective degree of freedom (edf), which measures the amount of smoothing, is large for all covariates. An edf of 1 indicates a straight line; thus, a model with edf = 1 is equivalent to a parametric linear model. A high edf (8-10) means the curve is highly nonlinear (Zuur et al. 2009). The evaluation of edf from Table 3 shows that a strong nonlinear relationship between the covariates and the stem radial measure. The inclusion of a smoothing term (nonparametric part) to the model is justifiable.

The effect of the covariates (Table 3) on stem radius might depend on the clone and season. Therefore, models that include clone, season and the interaction between clone and season in the parametric part were fitted for each covariate. Additive mixed models (AMM) with four smoothing terms for each of the climatic variables (one for each season) were fitted by including the interaction between clone and season in the parametric part of the additive mixed model. The likelihood ratio test comparing the model with an interaction and the model without an interaction term in the parametric part of these models were not significant (p ＞0.05). Therefore, models without the interaction effect of clone and season on the parametric part of each of the additive mixed models were used to study the effects of each covariate on the stem radius. The results for the effects of temperature, rainfall, relative humidity and wind speed showed a nonlinear relationship between stem radius and each of these covariates in autumn and winter. For summer and spring, the relationship can be explained by linear function. For solar radiation, there isnonlinear relationship with stem radius in autumn, and for the rest of the seasons, the relationship is linear. The type of relationship between stem radius and tree age also depends on season. The relationship is nonlinear in autumn and winter, and linear in summer and spring.

Table 1 Descriptive measures for stem radius

Table 2 Descriptive measures for climatic variables

Table 3 Comparison of models with linear trend and models with smooth terms

In the previous sections (Table 3), we applied a series of additive mixed models to stem radius for various covariates separately. In addition, tests for the presence of interaction between clone and season were made for different models.It was observed that the smoothers for tree age and each climatic variable also depend on season. The analysis made so far may help to see the effect of individual covariates on stem radius. To see the effect of more than one covariate on stem radius, it is essential to fit models that involve the smoothers for two or more covariates. This demands the application of model selection procedures. It is known that model selection with mixed models is complicated by the presence of fixed effects and random effects. The fixed effect structure and the random effect structure are dependent on each other and the selection of one affects the other. There are two strategies that are commonly used in a model selection process. These are the top-down strategy(Diggle et al. 2002) and the step-up strategy (West et al.2006). In the step-up strategy one starts with a limited model (e.g., few fixed and random effects) and then additional fixed effects and random effects are added based on statistical tests. In the top-down procedure, the initial model has one random intercept but with a model where the fixed component contains all explanatory variables and as many interactions as possible. This is called the beyond optimal model. Using the beyond optimal model, one can find the optimal component of the random effect (Zuur et al. 2009). The beyond optimal model is sometimes unrealistic due to many explanatory variables, interactions or numerical problems. In this paper, we followed the stepup approach.

Table 4 Parameter estimates of the additive mixed model with four smoothers of age and relative humidity and solar radiation (one per season in each case) with the interaction between season and clone included in parametric part(Maximum likelihood estimates)

For the step-up approach, both tree age and temperature were used to determine their effect on stem radius. The smoothed temperature for all four seasons was not significant. When the effect of the smoothed tree age was apparent in the model, the smoothed effect of temperature was not significant. An attempt to include temperature with one smoother for all four seasons in the model also shows that the smoothed temperature is not significant (edf = 1,p = 0.76). Instead of using temperature as a smoothed component of the AMM, we then used temperature in the parametric part of the AMM. A likelihood ratio test statistic was applied by including the interaction of temperature with season and the interaction of temperature with clone in the parametric part of the additive mixed model. In both cases, the interaction effect of temperature was not significant (p = 0.8 and 0.9 for the interaction with clone and season, respectively). When a likelihood ratio test was applied by including temperature in the parametric part of the additive mixed model, temperature was not important in explaining stem radial growth in the presence of the smoother for tree age in the model (p = 0.75). When the AMM that used the smoothers of tree age and any one of the climatic variables was fitted, the results showed significant smoothers for age in all cases. The smoothers for temperature, rainfall and wind speed did not appear to be significant. However, the smoothers for relative humidity for winter (p = 0.047) and the smoothers for solar radiation for winter (p ＜0.00001) and autumn(p = 0.0006) were significant. When an additive mixed model that included the smoothers of tree age, wind speed and solar radiation was fitted, all smoothers appeared to have an estimated edf of 1 (Table 4). The smoothers of tree age for all seasons (summer, autumn, winter and spring)were significant with small respective p values(＜0.00001). The smoothers for solar radiation were significant in autumn and winter with respective p values = 0.00171 and ＜0.00001. In all models used so far,the random intercept for each tree was used in combination with the assumption that residuals are normally distributed with mean 0 and constant variance. To validate the last model (that includes smoothers of tree age, solar radiation and relative humidity), the model validation graphs were plotted (Fig. 1). The lower right panel of the graphs shows a strong relationship between fitted and observed values of stem radius. The upper right panel shows that the assumption of constant variance was violated. The upper left and lower left panels show that there is some deviation from normality. In the plots of normalized residuals against covariates (tree age, solar radiation and relative humidity,clone and season) as part of the model validation process,there was no clear pattern as to the dependence of residuals on any of the covariates of tree age, solar radiation and relative humidity (Fig. 2). However, the spread of residuals depended on tree age. The spread of residuals also clearly depended on clone and season (Fig. 3), indicating that the variation in the data differed between seasons and clone. In addition, there was more variation in autumn and winter than in summer and spring, which violates the homogeneity assumption. Therefore, the assumption that the residuals are normally distributed with mean zero and constant variance was relaxed. A random slope was also used instead of random intercept.

Fig. 1 Basic model checking plots for the additive model with the smoothers of tree age solar radiation and relative humidity

All models were fitted using the REML estimation procedure and model comparison is made using the Akaike information criterion (AIC).

The model with random slope and different residual variance for each combination of clone and season had the smallest AIC and BIC (Table 5). The validation graph did not show any variation between seasons or between clones(Fig. 4). However, the plot of normalized residuals versus the fitted values (not shown) showed that there was still a certain degree of heterogeneity in the residuals. The last aspect of the modelling process was to allow for spatial or temporal correlation in the residuals. However, the attempt was not successful due to the complexity of the model.

Conclusion

Major climatic variables, nonclimatic variable tree age,tree clone, and the season were included in the study.Parametric techniques such as simple and multiple regression assume the linearity of the relationship between the response and the independent variables. However, the assumption of linearity is not realistic. Consequently, a semi-parametric approach was employed to explore the type of relationship between stem radius and the various covariates. Based on the fitted AMM for each individual covariate, we identified the type of relationship between stem radius and the covariates; stem radius has a nonlinear relationship with the covariates tree age and climatic variables (see Table 3, the significance of the smoothed terms).

Fig. 2 Plot of residuals versus tree age, solar radiation and relative humidity

Fig. 3 Normalized residuals plotted versus clone and season for the model that considers the smoothers of tree age, solar radiation and relative humidity

The effect of each climatic variable together with their interaction with clone was considered in the additive mixed model fitted to the data. To see the effect of more than one covariate at a time, we used a model selection procedure. Model selection for mixed models is much more complicated than for linear regression. For mixedmodels, the selection of the covariance structure is not easy due to computational and boundary value problems.A step-up model selection strategy (West et al. 2006) was applied to find a reasonable model that fits to the data,mainly because this strategy starts with the simplest possible model and is built up by including more covariates in the model and thus does not have much numerical problem. The model with random slope and different residual variance for each combination of clone and season had the smallest AIC and BIC; hence. It was the preferred additive mixed model for the data (Table 5). The results show that the GU clone has greater productivity potential than the GC clone as reported elsewhere (Galloway 2003; Melesse and Zewotir 2013a, b, 2015, 2017a, b, c). The effect of each covariate on stem radius was found to depend on season. A nonlinear relationship between stem radius and each covariate (tree age, temperature, rainfall, relative humidity and wind speed) was found in autumn and winter. For summer and spring, the relationship was captured by a linear function in a fully parametric model. The additive mixed models used in this paper are more informative than parametric models in the sense that they can identify the type of relationship (linear or nonlinear) in the data. This study used data from juvenile Eucalyptus trees,but the same technique will be used with more mature trees in the future.

Table 5 Akaike information criterion (AIC) and Bayesian information criterion (BIC) for models with different variance and correlation structure

Fig. 4 Normalized residuals plotted versus clone and season for the model with random slope and residual variance that varies with clone and season combination

AcknowledgementsThe authors are grateful to Dr. Valerie Grzekowiak and Dr. Nicky Jones for important comments and suggestions.