Linear mixed models with two longitudinal factors in the study of sugarcane dry root mass

There are agronomic experiments where measurements of a response variable are carried out over more than one longitudinal factor, for example, at different depths over time. These observations, made systematically in each experimental unit, can be correlated and might have heterogeneous variances at the different levels of the longitudinal factor. It was possible to model this correlation between repeated measures and the heterogeneity of variances by using mixed models. Thus, it was necessary to adapt some covariance structures that are common in experiments with only one longitudinal factor. The objective of this study was to use the class of linear mixed models to study sugarcane root dry mass. The experiment was the randomized complete blocks design and the parcels received four nitrogen doses. Repeated measurements were made over two longitudinal factors, one being qualitative ordinal (depths) and one being quantitative (distances from the planting line). It was possible to select a parsimonious covariance structure and another one to explain the average behavior of the responses through likelihood ratio tests, Wald tests, and using the AIC and BIC information criteria. The adjustment of the selected model was verified by using residual diagnostics graphs.


Introduction
In several research areas, it is common to conduct experiments in which repeated measurements of one or more variable responses are taken on successive occasions in the same experimental unit, which characterizes longitudinal experiments.Although experiments involving only one longitudinal factor are more common, there are cases where observations are made by considering two or more of these factors.In these studies, each experimental unit is associated with a vector of observations, called individual response profile, whose components are the observed values of the response variable at the levels of the longitudinal factor.
Considering that these repeated measurements are made systematically at each experimental unit, it is assumed that the observations made at different levels of the longitudinal factor are correlated and that the variances at different levels are heterogeneous (Lima, 1996).These characteristics of longitudinal data, which are not typical of usual regression models, can be included in the analysis model with the specification of linear (or non-linear) mixed models (Pinheiro;Bates, 2000;West;Welch;Galecki, 2007).Laird and Ware (1982), making use of the ideas of Harville (1977), proposed an approach that allows the use of different variances and covariance structures between and within experimental units involving the specification of two-stage linear mixed models when individual and population characteristics are identified.
In order to fit linear mixed models to a dataset involving a qualitative ordinal and a quantitative longitudinal factors, there is a need to adapt some variance and covariance structures that are common in studies with only one longitudinal factor.
The purpose of this work arose from the need to study a dataset by Otto (2007) in which nitrogen (N-urea) doses were applied in the furrow of a sugarcane plantation and the dry root mass (g dm -3 ) was observed at different depths (ordinal qualitative factor) and different distances from the sugarcane planting line (quantitative factor).
The objective of this work was to use the mixed linear model class to study the effects of the treatment factor in a randomized complete block design with repeated measurements along with two longitudinal factors, where one is qualitative ordinal and one is quantitative.

Matherials and Methods
The data analyzed are from Otto (2007) and come from experiments with sugarcane carried out at the Santa Adélia Plant, located in Jaboticabal, SP.Each experimental unit consisted of 48 furrows spaced 1.50 m apart (48 ×1.50 m = 72 m wide), with 15 m in length, totaling 1.080 m 2 per plot.The experimental design was a randomized complete blocks design with 4 treatments, referring to doses of 0, 40, 80 and 120 kg ha -1 N applied in the form of urea in the planting furrow, with 4 replications.The observed response variable was the sugarcane dry root mass (DRM), in g dm -3 .In order to collect the soil samples, a probe with internal diameter of 5.5 cm was used.The dry root mass was obtained after the samples were sieved (2 mm) and dried in an oven.The observations were made in July 2006 before the sugarcane harvest.DRM was evaluated over the planting row and at 30 and 60 cm from it, both to the right and left (Local: 0, 30 and 60) and also at two depths (Depth: 0 to 30 cm and 30 to 60 cm), according to the Figure 1.It can be seen in Figure 2 a greater variability in sugarcane dry root mass at Depth 0-30 cm and in the planting row (Local = 0).It is also noted that there is a decrease in DRM as the distance from the planting row increases.Regarding the variability, a decrease is noted as the distance from the same row increases.Thus, the closer it is to the sugarcane, the higher concentration of DRM can be found, however, not all observations made close to the sugarcane showed this high concentration.
To fit a mixed model to these data, the approach of Laird and Ware (1982) was used, which uses the linear model, where, y ij (b × 1) is the response profile of the (ij)-th individual, being block i, treatment j, for i = 1, . . ., r and j = 1, . . ., a; X ij (b × p) is a known and specification rank matrix p < b associated with the vector β(p × 1) of unknown subpopulational parameters, where b is the number of longitudinal factor levels and p the number of fixed effect parameters; Z ij (b × q) is a known and specification matrix, with a complete column rank and associated with the vector of random effects u ij (q × 1) of individual differences around the population values, and ε ij (b × 1) is a vector of random errors.
The model (1) can be formulated in two stages, and in the first stage, for each experimental unit (ij), there is: where, R ij is known as conditional dispersion matrix and is associated with the conditional error In the second stage, it is assumed that u ij ∼ N(0, G) is independent of ε ij , thus it is obtained the marginal (or unconditional) model: where, the matrix is called marginal dispersion matrix and is associated with the marginal error e ij = y ij -X ij fi, which is given by the sum To establish the fixed and random effects to compose the mixed model, it was used the top-down strategy (Verbeke;Molenberghs, 2000), which began by fitting the maximal model with all possible fixed effects terms of interest in order to explain thoroughly the systematic variation in the responses.Subsequently, we sought a structure for the random effects (covariance matrix G) to be included in the model.The remainder of the variability, which was not explained by the variation betweenindividuals, is due to the within-individual experimental error, selecting an appropriate covariance structure for the covariance matrix R ij .In comparisons involving the fixed effect parameters, the likelihood ratio test for nested models (Verbeke; Molenberghs, 2000; West; Welch; Galecki, 2007) and the Akaike (Akaike, 1974) and Bayesian (Schwarz, 1978) information criteria were used.After visual analysis of individual DRM profile plots along the quantitative longitudinal factor levels, for each qualitative longitudinal factor level (Figure 2), and still observing the average profiles by dose and by depth of Figure 3, the maximal model ( 4) was defined for the fixed part involving a second degree polynomial for each of the nitrogen doses: for i = 1, 2, 3, 4; j = 1, 2, 3, 4; k = 1, 2 and l = 1, 2, 3, where, E(y ijkl ) is the average DRM of the parcel, in the ith Block, that received the jth Nitrogen dose, at the kth Depth and lth Location; B i is the fixed effect of the ith block; β 0jk , β 1jk and β 2jk are, respectively, the intercept, the first degree coefficient and the second degree coefficient of the parabola associated with the parcels (jk).In this notation:

for all combinations (jk).
The marginal dispersion matrix (V ij ) defined in (3) can be modeled by specifying random effects (Z ij ), with its corresponding covariance structure between individuals (G) and the specification of a structure for within-individual measurements (R ij ).Firstly, it is aimed to model the (positive) correlation between experimental units and the heterogeneity of variances.In order to achieve a better fit to the structure of sample variances and correlations, the process is completed by verifying the need for a more complex structure for the matrix R ij .The process of choosing the random effects and the respective covariance structure started with the inclusion of a random effect in all the coefficients of the polynomial curve and a completely parameterized covariance structure.
In Table 1 some interesting structures for the within-individual covariance matrix (R ij ) are presented, they are adapted for the case of two longitudinal factors.These structures were defined by assuming that the three levels of Locations (second longitudinal factor) are nested within each of the two levels of Depths (first longitudinal factor).Such structures were created by using the pdMat, varFunc, and corStruct classes arguments of the lme function of the nlme library (Pinheiro;Bates, 2000) of the R program (R Core Team, 2022).Other structures can be created by combining the arguments of these classes, depending on the variability and correlation characteristics of the data.

Estimation and selection methods
In the process of choosing the random part of the mixed model, the restricted maximum likelihood (REML) method was used to estimate the parameters of the models, since the fixed part of the model was kept unchanged.Models with nested random effects structures were compared by using the likelihood ratio test (LRT) and the Akaike (AIC) and Bayesian (BIC) information criteria.According to West, Welch and Galecki (2007), for cases of comparisons in which the parameters lie on the boundary of the parameter space, the LRT statistic follows a mixture of χ 2 distributions.For example: when the inclusion of a random effect in the model is tested, the null hypothesis (H 0 ) imposes that the variance of this random effect must be null.Negative values for the test statistic are not allowed.In this case, the descriptive level of the test was calculated by where χ 2 H 0 and χ 2 H 1 are the quantiles of the chi-squared distribution with degrees of freedom given by the random effects numbers of the models under the assumptions H 0 and H 1 , respectively, and LRT is the likelihood ratio calculated value.
After choosing the random part of the mixed model, the significance of the parameters of the fixed part was tested, we fit the models by the maximum likelihood (ML) method, by using the Wald, LRT tests if the models are nested, and the AIC and BIC information criteria if models are not nested.The fixed effects selection process was ended when a parsimonious model was obtained, which managed to explain with a smaller number of parameters all the important characteristics of the phenomenon being studied.
The diagnosis of the fitted models was performed by using a plot of standardized conditional residuals versus predicted values in order to verify homoscedasticity of these residuals; plot of standardized conditional residuals versus observations' indexes to detect outliers; the normality assumption was verified by plot of norm quantiles versus standardized conditional residuals and by histogram of these residuals (Singer; Nobre; Rocha, 2018), furthermore, it was confirmed by nonparametric Shapiro-Wilk test (Shapiro;Wilk, 1965).
Finally, to visualize the goodness of fit of the final selected model, scatterplots of the observed data were constructed, with the mean curves fitted and the individual curves predicted.

Results and Discussion
While evaluating the descriptive statistics of dry root mass (Table 2) for the main effects of treatment and longitudinal factors it was observed that the mean and standard deviation of the DRM in the doses 0 and 80 kg ha -1 of N-urea were very close and lower than the responses of the doses 40 and 120 kg ha -1 , and also, that the mean and standard deviation of the DRM decrease with increasing depth or distance from the planting row (local).
Compound symmetry with heterogeneity: Equal variances at different locations at the same depth, different variances at the two depths, equal covariances between depths at different locations, and equal correlation between depths.Number of parameters: 4.
Variance components with heterogeneity: Different variances at the two depths, but equal in the three locations evaluated at the same depth.Number of parameters: 2.

Random effects -matrix G
In addition to the intercept, the inclusion of random effects to the other parameters of the curves was tested, as indicated in Table 3.
In this initial phase, the structure R = Iσ 2 was used for the within-individual covariance matrix.Models were fitted by using the REML method and compared by using the likelihood ratio test (LRT).The results (Table 4) were favorable to the m3 model, which proposes random effects for the intercepts of the two curves (Depths 0-30 cm and 30-60 cm) and for the coefficients of degree 1 and degree 2 of the curve relative to Depth 0-30 cm.The analysis did not indicate the need to include random effects to the coefficients of the first and second degrees of the curve relative to Depth 30-60 cm, due to the low variability of the responses in different locations.The structure of the random effects matrix was where, D1 -Depth 0-30 cm and D2 -Depth 30-60 cm.

Covariance matrix (R)
The suggested structures are in Table 5, where, • m5, significant difference for Depth level; • m6, significant interaction of Depth and Location levels; • m7, homogeneity for levels 30 and 60 for Locations; • m8, Dose levels differ for Depth levels and Location levels; • m9, Doses 0 and 80 are equal and the others are different.
equal for doses 0 and 80 and different for other levels The results of the comparisons between the models m5, m6, m7, m8 and m9, described in Table 5, with the model m3 chosen in the previous section, are in Table 6.In the comparisons contained in Table 6, it was observed that among the nested structures (m3, m5 and m6) the LRT indicated the choice of the structure m6, which proposes different variances at each depth and at each location.The comparison of the m6 and m7 models, which are not nested, was favorable to the m7 model, which presented the lowest AIC and BIC values.Comparing the m7 and m8 structures by the LRT, the m7 structure was favorable, and it allowed to conclude that the DRM data for the different doses have equal variances.And lastly, in the comparison of the m7 and m9 structures by the AIC and BIC values, the m7 model was favorable, which confirms that the variances of the DRM in the different doses were considered equal.
Other simpler (CS, AR(1)) or more complex (UN) intra-individual structures were compared to the diagonal structure with heterogeneity (diag w/ het) of the m7 model (Table 7), but all comparisons were favorable to the m7 structure, which has a smaller number of parameters.Therefore, the selected covariance structure was G m3 (6) and R m7 (Table 5).

Fixed effects
To verify the equality of the regression coefficients of the fixed maximal model, and the chosen covariance structure, the Wald test was used, with the model fitted by the maximum likelihood (ML) method.From the results presented in Table 8, it was possible to conclude that the regression coefficients of degrees 0, 1 and 2 of the maximal model ( 4) are significantly (p < 0.05) non-null.In Table 9 are the structures that were compared in the sequence of the analysis.These comparisons are shown in Table 10.The m10 structure which fits a curve for each of the four Dose levels for the Depth 0-30 , and only one curve for the Depth 30-60 was compared with the m7 structure, which considers a curve for each Dose level and for each Depth level.According to the information criteria, the coefficients of degrees 0, 1, and 2 for Depth 30-60 do not differ significantly from each other.
Next, the structure m11 was fitted, considering only one curve for the Depth 0-30 .Comparing m10 with m11, the equality of the coefficients of degrees 0, 1 and 2 of the four curves for the Depth 0-30 should not be considered, since the comparison by the AIC indicates that the model m10 fits better than m11.On the other hand, the comparison by BIC, which penalizes the model with the highest number of parameters, shows that the best fit is the m11.Therefore, in order to compare the equalities of the means of combinations of the Dose levels, it was considered that the structure m11 was rejected by the AIC.
Thus, the structure m10 was compared to the structure m12, in which we consider the equality of the coefficients of the curves referring to Doses 0 and 80 kg ha -1 N for the Depth 0-30 .In this comparison, the equality of these two curves must be considered, since the AIC and BIC criteria indicated that the structure m12 has a better fit in relation to the m10.
Thus, the structure m12 was compared to m13 in which, in addition to the equality of the curves referring to Doses 0 and 80, in the last one it is also considered the equality of the curves referring to Doses 40 and 120 kg ha -1 N to Depth 0-30 .For this comparison, according to the AIC and BIC information criteria, there are only two curves for the Depth 0-30 and only one for the Depth 30-60 .
Furthermore, the significance of the quadratic coefficient for the Depth 30-60 (m13) in relation to a straight line (m14) was tested.In this comparison, it resulted that the quadratic coefficient is significant to the model, therefore, the m13 structure is the one that best fits to the DRM data.

Diagnostic of the selected model -m13
Figure 4 shows the residuals diagnostic plots.In 4(a), it is observed that there is a greater variation in the amplitude of the standardized conditional residuals for the predicted values closer to zero.However, since there is no point outside the -3 to 3 range.In 4(b), it can be observed that there are no outliers.
Regarding the normality of the residuals, it can be seen in the plots 4(c) and 4(d) that the standardized conditional residuals of the selected model (m13) are very close to the normal distribution, which is confirmed by the p-value = 0.221 of the Shapiro-Wilk test.

Final Model -m13
Parameter estimates for this model are presented in (7), the random effects matrix G in (8), and the within-individual covariances matrix (R ij ) in (9).R ij = diag 0.012, 0.002, 0.002, 0.022, 0.004, 0.004 (9) Figure 5 shows the observed values of sugarcane DRM, as well as the average and individual curves for the final model (m13).In this same figure, it can be observed that at Depth 0-30 cm the Nitrogen doses that result in a higher sugarcane DRM are 40 and 120 kg ha -1 N. The curve fitted for the Depth 30-60 cm was below the two curves for the Depth 0-30 cm, which indicates that regardless of the nitrogen dose applied, the greater the depth, the lower the sugarcane DRM always is.It can also be observed that in the planting row the DRM decreases more sharply at the Depth 0-30 cm, whereas when it is beyond 30 cm from the planting row, the DRM at both depths and with four nitrogen doses are very similar.
It can also be noted in Figure 5 that the selected covariance structure captured the variability of the data properly, which indicates that the final linear mixed model also explained well the behavior of individual responses and average responses as a function of increasing nitrogen doses, depths and locations.

Conclusions
The use of linear mixed models in the analysis of dry root mass data, from an experiment involving two longitudinal factors, was quite effective, as it managed to properly explain the behavior of average response and the complex structure of correlations between the levels of these factors and heterogeneous variances.
The functions of the nlme library from the R package were very important in the analysis, for they provided the adjustment of several covariance structures involving the two longitudinal factors, the comparison between them, and the construction of excellent quality graphics used in the diagnostics.

Figure 1 .Figure 2 .
Figure 1.Illustration of soil sample collection for plant cane

Figure 3 .
Figure 3. Average profiles of sugarcane dry root mass per dose at each depth

Figure 5 .
Figure 5. Observed DRM values, fitted mean curves and predicted individual curves in the final model -m13

Table 1 .
Structures for the matrix R ij adapted to two longitudinal factors, a qualitative factor with 2 levels (depth) and a quantitative factor with 3 levels (local)

Table 3 .
Initial structures for inclusion of random effects, considering the maximal structure (4) for fixed effects, for sugar-

Table 4 .
Initial comparison of models aiming to include random effects by LRT

Table 5 .
Within-individual covariance structures (R) for sugarcane DRM data

Table 6 .
Comparison of m3, m5, m6, m7, m8 and m9 models by LRT and by the AIC and BIC information criteria

Table 8 .
Wald test for fixed-effect parameters of the maximal model (m7)