INTRODUCTION
In tropical livestock systems, poor grassland productivity is one of the most important limitations due to the adaptability and persistence in these environments. However, new cultivars of the genus Brachiaria have been released to the market as options to overcome the problems observed in traditional forages, thus providing better fodder options^{1}. Brachiaria brizantha is one of the superior introducing grasses that has been adapted and known by farmers. These grasses are compatible with the tropical climate and are tolerant of various types of soil, including acidic soils^{2}. Knowledge of the distribution of production and quality of forage during the year is a tool to plan utilization^{3}. Forage production in the tropics is related to the distribution of rainfall, which is generally asymmetrically distributed. This fact raises the need to use more complex models, especially when multiple comparisons are made and there are very large deviations from normality, the risk of error increases^{4,5}.
Generalized Additive Models for Location, Scale and Shape (GAMLSS) are semiparametric regression type models. They are parametric, in that they require a parametric distribution assumption for the response variable and “semi” in the sense that the modeling of the parameters of the distribution, as functions of explanatory variables, may involve using nonparametric smoothing functions. GAMLSS were introduced by researchers^{6} as a way of overcoming some of the limitations associated with the popular generalized linear models, GLM and generalized additive models, GAM^{7,8}. In GAMLSS the exponential family distribution assumption for the response variable (y) is relaxed and replaced by a general distribution family, including highly skew and/or kurtotic continuous and discrete distributions^{9}. The systematic part of the model is expanded to allow modeling not only of the mean (or location) but other parameters of the distribution of y as, linear and/or nonlinear, parametric and/or additive nonparametric functions of explanatory variables and/or random effects^{6}. Hence GAMLSS is especially suited to modeling a response variable that does not follow an exponential family distribution, (e.g., leptokurtic or platykurtic and/or positive or negative skew response data, or overdispersed counts) or which exhibit heterogeneity, (e.g., where the scale or shape of the distribution of the response variable changes with explanatory variables(s))^{9}.
In this paper, a GALMSS model for skewed data, where the exponential family assumption is relaxed and replaced by a very general distribution family is considered. Within this new framework, the systematic part of the model is expanded to allow not only the mean (or location) but all the parameters of the conditional distribution of y to be modeled as parametric and/or additive nonparametric (smooth) functions of explanatory variables and/or randomeffects terms. The objective of this research is to fit a GAMLSS model on accumulated dry matter data from Brachiaria brizantha using a model selection algorithm.
MATERIALS AND METHODS
Study area: The study was carried out at a farm in Portuguesa state, Venezuela (Fig. 1) from November, 2016October, 2017). The farm is located between 10080201004802 West Longitude and 425400427420 North Latitude. The response variable ‘y’ is the accumulated dry matter (kg ha^{–}^{1}). The explanatory variable ‘x’ is the interval between cutting (21, 28, 35 and 42 days).
BoxCox power exponential distribution: The BoxCox power exponential distribution (BCPE) is defined by (for details see^{11}).
Let Y be a positive random variable having a BoxCox power exponential distribution, denoted by , defined through the transformed random variable Z given by:
for 0<Y<∞ where y and where the random variable Z is assumed to follow a standard power exponential distribution with power parameter, τ>0, treated as a continuous parameter (Parameterization (1) assumes a standard normal distribution for Z)^{11}.
The probability density function of Z, a standard power exponential variable, is given by:
for ∞<z<∞ y τ>0, where c^{2} = 2^{–}^{2/}^{τ}Γ(1/τ)[Γ(3/τ)]^{–}^{1}. This parameterization ensures that Z has mean 0 and standard deviation 1 for all τ>0. Note that τ = 1 and 2 correspond to the Laplace (i.e., twosided exponential) and normal distributions respectively, while the uniform distribution is the limiting distribution as τ→∞^{12}.

Fig. 1:  Relative location of the Portuguese state, Venezuela 

^{Source: Seijas10} 
From 2, the probability density function of Y, a BCPE (μ, σ, υ, τ), random variable, is given by:
Modeled data using GAMLSS: Let y^{T} = (y_{1}, y_{2},..., y_{n}) be the vector of the response variable observations. Also, for k = 1,2,..., p, let g_{k}(.) be a known monotonic link function relating θ_{k} to explanatory variables and random effects through an additive model given by:
where, θ_{k} and η_{k} are vectors of length n, e.g., θ^{T}_{k} = (θ_{1k}, θ_{2k},..., θ_{nk}), β^{T}_{k} = (β_{1k}, β_{2k},..., β_{j’kk}) is a parameter vector of length J’_{k}, X_{k} is a known design matrix of order n×J’_{k}, Z_{jk} is a fixed known n×q_{jk} design matrix and γ_{jk} is a q_{jk}dimensional random variable. Model (4) is called the GAMLSS.
The vectors γ_{jk} for j = 1,2,..., j_{k} could be combined into a single vector with a single design matrix Z_{k}.
If, for k = 1,2,..., p, J_{k} = 0 then model (4) reduces to a fully parametric model given by:
If Z_{jk} = I_{n}, where I_{n} is an n×n identity matrix and γ_{jk} = h_{jk} = h_{jk} (x_{jk}) for all combinations of j and k in a model (4), this gives:
where, x_{jk} for j = 1,2,..., J_{k} and k = 1,2,..., p are vectors of length n. The function h_{jk} is an unknown function of the explanatory variable x_{jk} and h_{jk} = h_{jk} (x_{jk}) is the vector that evaluates the function h_{jk} at x_{jk}. The explanatory vectors are assumed to be known. We call the model in equation (5) the semiparametric GAMLSS. Model (6) is an important special case of model (4). If Z_{jk} = I_{n} and γ_{jk} = h_{jk} = h_{jk} (x_{jk}) for specific combinations of j and k in a model (4), then the resulting model contains parametric, nonparametric and randomeffects terms.
The first two population parameters θ_{1} and θ_{2} in a model (4) is usually characterized as location and scale parameters, denoted here by μ and σ, whereas the remaining parameter(s), if any, are characterized as shape parameters, although the model may be applied more generally to the parameters of any population distribution.
For many families of population distributions a maximum of two shape parameters ν (=) and τ (=θ_{3}) suffice, giving the model:
The GAMLSS model (4) is more general in that the distribution of the dependent variable is not limited to the exponential family and all parameters (not just the mean) are modeled in terms of both fixed and random effects, for details see Rigby and Stasinopoulos^{11}.
Model for the four parameters of the BoxCox power exponential distribution: The parameters (μ, σ, ν, τ) of the BoxCox power exponential distribution may be modeled as functions of many explanatory variables using the generalized additive model for location, scale and shape (denoted by GAMLSS) of Rigby and Stasinopoulos^{11}. Here we consider a special case of the GAMLSS model where there is a single explanatory variate X. Given X = x; Y is modeled by a BoxCox power exponential variable, BCPE (μ, σ, ν, τ), with probability density function f_{Y}(y), defined by (4) where the parameters μ, σ, ν and τ are modeled as smooth nonparametric functions of x, i.e. Y∼BCPE (μ, σ, ν, τ) where:
and for k = 1, 2, 3, 4, g_{k}(.) are known monotonic link functions, usually the identity for μ and ν and the log for σ and τ and h_{k} (x) are smooth nonparametric functions of x.
For i = 1,2,..., n, given X = x_{i}, observations Y_{i} are assumed to be independent BCPE (μ_{i}, σ_{i}, ν_{i}, τ_{i}) variables with probability density functions f_{Yi}(y_{i}) obtained from (3) and parameters obtained from (8). This model is appropriate for independent observations of Y (e.g. crosssectional data) rather than correlated observations (e.g. longitudinal data), for details see Rigby and Stasinopoulos^{11}.
Model estimation and selection: The nonparametric functions h_{k} for k = 1, 2, 3, 4 are estimated by maximizing the penalized loglikelihood function l_{p} defined by:
where, h”_{k}(u) is the second derivative of h_{k}(u) to u and is the loglikelihood function of the data and l_{i} is the loglikelihood function of observation y_{i} from a BoxCox power exponential distribution, BCPE (μ_{i}, σ_{i}, ν_{i}, τ_{i}), obtained from (4).
The penalized loglikelihood function (9) is maximized iteratively using either the RS or CG algorithm of Rigby and Stasinopoulos^{11}, which in turn uses a backfitting algorithm to perform each step of the Fisher scoring procedure, requiring the loglikelihood of the data, I_{d} and first and expected second derivatives to μ, σ, ν and τ.
A general criterion for model selection is the generalized Akaike Information Criterion (GAIC), obtained by adding to the fitted deviance a fixed penalty # for each effective degree of freedom used in a model, i.e. GAIC (#) = where df denotes the total effective degrees of freedom used in the model, denotes the fitted deviance and denotes the fitted loglikelihood.
The total effective degrees of freedom df combines the effective degrees of freedom used in the smooth functions h_{1} (x) to h_{4} (x) in (8), for modeling μ, σ, ν and τ, denoted by df_{μ}, df_{σ}, df_{ν} and df_{τ}, respectively. Each effective degrees of freedom (e.g. df_{μ}) is defined by the trace of the corresponding smoothing matrix in the fitting algorithm, which is in turn directly related to the corresponding smoothing parameter (e.g. λ_{1}), see Rigby and Stasinopoulos^{11} for more details. The model with the smallest value of the criterion GAIC (#) is then selected. The Akaike information criterion (AIC) and the Schwartz Bayesian criterion (SBC) are special cases of the GAIC (#) criterion corresponding to # = 2 and # = log (n), respectively. Let BCPE (df_{μ}, df_{σ}, df_{ν}, df_{τ}, λ) represent the BCPE model (8), where the first four values inside the brackets denote the total effective degrees of freedom used in the smooth nonparametric functions h_{1} (x) to h_{4} (x) for modeling μ, σ, ν and τ, respectively and the fifth value λ denotes the power transform parameter in the transform x = x var^{λ}, where xvar is the explanatory variate recorded in the data set.
For cumulated dry matter data set the threestep procedure for selecting a model of form (8) proposed by Rigby and Stasinopoulos^{11} is used:
Initial choices:
• 
Choose link functions g_{k} (.) For k = 1, 2, 3, 4 in (4) 
• 
Choose an initial parameter value λ_{0} for λ in the transformation x = x var^{λ}, assuming x var>0 for all observations. The initial value λ_{0} is chosen to provide an approximate constant absolute gradient between y and x. (A shifted (or offset) power transformation of x var, i.e., x = (x var+δ)^{λ} could be considered if x var takes negative values) 
• 
Choose a single penalty # for each effective degree of freedom used in the models for μ, σ, ν and τ 
Model selection:
• 
From model BCPE (1, 1, 1, 1; λ_{0}) (where the first four values df_{μ} = df_{σ} = df_{ν} = df_{τ} = 1 indicate fitting a constant for parameters ), apply forward selection of to give the chosen value which minimizes criterion GAIC (#) 
• 
In model BCPE (df_{μ}_{1}, 1, 1, 1, λ), estimate λ, (e.g. using a grid search over values of λ) to give the value λ_{1} which minimizes GAIC (#) 
• 
From model BCPE (1, 1, 1, 1; λ_{1}), apply forward selection of df_{σ} to give chosen value df_{σ1} which minimizes GAIC (#) 
• 
From model BCPE (df_{μ}_{1}, df_{σ1}, 1, 1, λ_{1}) create a twoway table of values of GAIC (#) for combinations (df_{ν}, df_{τ} ) and search for the combination (df_{ν1}, df_{τ1}) which minimizes GAIC (#) 
Finetuning of model BCPE (df_{μ}_{1}, df_{σ1}, df_{ν1}, df_{τ1} λ_{1}):
• 
Finetuning of df_{σ}, by changing df_{σ1}, in steps of 1, if the value of GAIC (#) decreases 
• 
Similar _ne tuning of df_{μ} 
• 
Finetuning of λ by reestimating λ as in step 2 (b) 
The R code to use the algorithm described above is shown in the Appendix.
RESULTS AND DISCUSSION
Dry matter of Brachiaria brizantha: In Fig. 2a, the distribution of the dependent variable (accumulated dry matter) of B. brizantha showed asymmetry and kurtosis. Figure 2b the box plots show an asymmetric distribution of the accumulated dry matter data for intervals between cuts of 35 and 42 days. Figure 2c the densities associated with the accumulated dry matter data of B. brizantha show leptokurtic distributions for the 4 intervals between cuts. The data of Table 1 shows how the accumulated dry matter of B. brizantha increases as the interval between cuts increases. The data of
Appendix:  R code to fit a GAMLSS model on dry matter data in Brachiaria brizantha 


Fig. 2(ac): 
Dry matter distribution of Brachiaria brizantha 

^{(a) Histogram, (b) Boxplot and (c) Density at four intervals between cutting} 

Fig. 3:  BoxCox Power Exponential Distribution (BCPE) fitted on Brachiaria brizantha dry matter data at four intervals between cuttings 
Table 1:  Yield (kg) of dry matter of Brachiaria brizantha at four intervals between cuttings 

Table 2:  Assumptions on dry matter data from Brachiaria brizantha at four intervals between cutting 

Table 3:  GAMLSS model adjusted on dry matter data of Brachiaria brizantha at four intervals between cutting 

Table 2 shows the nonnormality or asymmetry (p<0.05) of the accumulated dry matter distribution and serial autocorrelation of errors (p<0.05). From Fig. 3 and Table 3, the chosen BCPE model provides a substantially better fit than the previous models both in terms of the GAIC (865).
Table 4:  BoxCox power exponential distribution (BCPE) fitted on Brachiaria brizantha dry matter data at four intervals between cutting 


Fig. 4(ac): 
Fitted values for the parameters (a) μ,(b) σ and (c) ν in the final chosen BCPE model on dry matter data in Brachiaria brizantha 
The fitted models for the μ, σ and in the chosen BCPE model have plotted in Fig. 4ac, respectively. The estimated parameters for BoxCox power exponential distribution (BCPE) fitted on Brachiaria brizantha dry matter data at four intervals between cuttings are shown in Table 4. The parameters μ, σ, ν and τ of the BCPE distribution may be interpreted as relating to the location (median), scale (coefficient of variation), skewness (power transformation to the symmetry) and kurtosis, respectively^{9}. The fitted value for τ in Equation (13) indicates leptokurtosis as τ>3. From Table 5 for the different intervals between cutting, the distribution of dry matter from B. brizantha is asymmetric as v_{i}>0. From Eq. 10 and Fig. 4a, it is observed that the longer the interval between cuts, the average yield of dry matter increases. Similar results are reported by Hare et al.^{13} in a study on the effect of cutting interval on yield and quality of three brachiaria hybrids in Thailand. From Fig. 4b it is observed that the longer the interval between cuts, the variability in the dry matter yield decreases and from Fig. 4c it is observed that for intervals between the cutting of 2128 days the asymmetry of the distribution of dry matter yield decreases, while for intervals between the cutting of 3542 days it increases.
This indicates that cutting stimulates regrowth of the plant and consequently its dry matter production since a significant increase in dry matter is observed when going from the first to the third cut. There is enough information in the literature to explain the effect of the cutting stimulus^{13} or grazing^{1,14} in the dry matter production of pastures. These variations in the average dry matter production between cuts of the genotypes as a whole were the product of the environmental effect (heterogeneous distribution of rainfall within the same dry season), rather than by effects of the genotype itself (p>0.05). In this regard, some studies show that environmental effects are often more important than the effects of the genotype on the productive traits of forage plants^{15,16}. In this way, the modeling of the accumulated dry matter using the BCPE distribution allows to show the effect of the interval between cut on the accumulated dry matter and to make decisions in the management of this species that allow better yields.
Table 5:  Estimated parameters from the GAMLSS model on the accumulated dry matter from B. brizantha at four intervals between cuttings 

In the parameterization of the BoxCox Power Exponential Distribution (BCPE) fitted on Brachiaria brizantha dry matter data at four intervals between cutting with parameters used in the GAMLSS methodology, its expectation is given by: E[Y] = μ and its variance by Var[Y] = σ^{2}μ^{2}:
CONCLUSION
GAMLSS allowed flexible modeling of both the distribution of the dry matter yield of B. brizantha and the dependence of all the parameters of the distribution on intervals between cutting. The BCPE distribution provides a flexible model for the dry matter yield of B. brizantha, exhibiting both skewness and leptokurtosis. For the dry matter yield from Brachiaria brizantha, exhibiting skewness and platykurtosis, the BCPE distribution, provided the best fit. Finally, the results of this research suggested that the interval between cutting has an effect that is reflected in the average yield of dry matter of B. brizantha, that is, the yield of dry matter increases as the interval between cutting increases. . Similarly, the interval between cuts affects the skewness and the kurtosis of the distribution.
SIGNIFICANCE STATEMENT
This study uncovers the potential uses of GAMLSS models in the dynamics of dry matter production that may be beneficial to pasture and forage researchers. This study uncovers the potential uses of GAMLSS models in the dynamics of dry matter production that may be beneficial to pasture and forage researchers. This study will help the pasture and forage researcher to uncover critical areas of accumulated dry matter production of Brachiaria brizantha that many researchers were unable to explore. Therefore, a new theory can be reached on modeling the dynamics of accumulated dry matter production of B. brizantha.