In vivo studies of the nutritional value of foods are time consuming. They require a sufficient number of animals and large amounts of plant material, and are, therefore, economically expensive (Posada and Noguera 2005). However, in vitro techniques are less expensive, performed in less time, require only small amounts of plant material as substrate, and allow greater control of experimental conditions (Makkar 2005). In this regard, Theodorou et al. (1994) proposed a method to study the fermentation kinetics of substrates from the gas production generated by their in vitro fermentation, which is estimated at different times from the gas pressure measured inside glass bottles. The volume is assumed constant for all bottles and over time.
In general, the models to describe in vitro gas production kinetics are non-linear models (NLM). They can be classified according to their structure in exponential, sigmoidal and multicompartmental or multiphase models (Muro et al. 2017). However, a frequent problem when modeling in vitro gas production data is that residues are correlated with each other. One of the possible causes is that, in order to obtain in vitro gas production, repeated measurements are carried out over time on the same experimental unit and a specific structure of variances is generated. For this reason, mixed NLMs constitute an alternative for analysis. With them, the correlation structure can be modeled, directly or through random variables. Other benefits of the mixed approach are that it does not require data distribution to be normal, help to control heterogeneity, improve the statistical adjustment criteria and have positive effects in the analysis of NLM, fitted to longitudinal data (Gómez and Agüero 2020).
Mixed models gain value in the agricultural and livestock field because they are applied to the analysis of longitudinal data, multi-environment trials and growth curves (Bandera and Pérez 2018). The use of linear models with a mixed approach is common in the analysis of variance of in vitro gas production experiments with repeated measures (Gómez et al. 2019). However, little is known about the use of mixed NLMs to describe in vitro gas production performance.
The main objective of this study was to evaluate the performance of mixed NLMs and NLMs in the description of in vitro gas production kinetics of ruminant feeds.
Materials and Methods
Experimental procedure. Data was selected from an experiment that was carried out at the Institute of Animal Science in 2019. In this study, the effect of including sweet potato tuber and Lactobacillus pentosus LB-31 on the nutritional value of the mixed silage of Cenchrus purpureus x Cenchrus glaucum (hybrid OM-22) + Moringa oleifera was evaluated (Rodríguez et al. 2019). The experiment used in vitro gas production (GP) technique, described by Theodorou et al. (1994). A completely randomized experimental design was applied, with factorial arrangement (3 x 2). The factors were three levels of inclusion of the tuber (0, 25 and 50 %) and the use or not of the bacterial strain. GP was measured at 3, 6, 9, 12, 15, 18, 21, 24, 29, 48, 72, 77 and 144 h and three repetitions were performed per treatment. Although the experiment analyzed six silage variants, for the present study, GP data from the fermentation of mixed silage + 25 % of sweet potato tuber + the microbial additive (L. pentosus LB-31) were used.
Statistical analysis. The sigmoidal equations shown in table 1 were modeled. Three analyzes were performed for the evaluation. In the first, only the random effect of residuals was considered. In the second, a random effect was added to parameter b, as it was the one that showed the greatest variability among the incubation bottles. This random effect is related to the individuality of each incubation bottle. With its inclusion in the model, an attempt was made to control the source of variability that each individual subject might have.
Analysis | Logistic model |
Gompertz model |
|
---|---|---|---|
1: GP (t) = | (b)/(1+ exp(2-4*c*(t-L)))+e | (b)*exp(-k*exp(-r*t))+e | (b)/(1+(T/t)F)+e |
2: GPi (t) = | (b+bi)/(1+exp(2-4*c*(t-L)))+e | (b+bi)*exp(-k*exp(-r*t))+e | (b+bi)/(1+(T/t)F)+e |
3: GPij (t) = | (b+bij)/(1+exp(2-4*c*(t-L)))+e | (b+bij)*exp(-k*exp(-r*t))+e | (b+bij)/(1+(T/t)F)+e |
GP: gas production at time t (mL g-1 incubated OM), c: GP rate (h-1), t: fermentation time (h), L: Lag phase (h), b: asymptote when t = ∞ (mL g-1 incubated OM), k: relative GP rate (h-1), r: constant factor of microbial efficiency (h-1), T: time in which half of the GP (h) is reached, F: constant that determines shape and characteristics of curve profile, i: incubation bottle, j: incubation time, e: error term, bi: asymptotic random effect associated with bottle “i”, bij: asymptotic random effect associated with bottle "i" at time "j".
In the previous analyzes, the fact that GP data are not taken with total randomness was neglected, but that measurements are made repeatedly on the same experimental unit or incubation bottle, violating the assumption of independence of observations. Therefore, in the third analysis, it was considered that there is a random effect associated with asymptotic GP, and that it varies according to the bottle and the bottle-time incubation time.
The estimation of parameters was carried out with the NLMIXED procedure of SAS 9.3 (2013). To apply this procedure, it was assumed that observations were independent, and that there was a joint probability density function, which links observations and vectors of random effects. Random effects were assumed to be normal, with zero mean and constant variance. They were estimated by the empirical Bayes method (Gaver and O'Muircheartaigh 1987).
The parameters obtained with the NLMIXED proc were substituted in the MODEL proc. In this way, White test (Pérez et al. 2020) was executed to assess the homogeneity of variances of residuals, and Durbin-Watson (DW) test was used for determining the presence of autocorrelation (Hurtado et al. 2021). With the help of the CAPABILITY proc, the Kolmogorov-Smirnov (De Roa 2020) normality test was applied. In addition, the SPSS 22 (2013) program applied the non-parametric run test for randomness to residual. This test allowed to contrast the hypothesis of a random order versus a trend alternative.
The selection of the model with the best fit was made based on the variance of the error (σ2), the coefficient of determination adjusted to the degrees of freedom (R2aj), the significance of the parameters, Akaike Information Criterion (AIC) and the Bayesian Information Criterion (BIC) (Tiago et al. 2017). Models with lower σ2, AIC and BIC values were considered to have a better fit, unlike R2aj, in which higher values are preferred. The fulfillment of hypotheses on the residual component was assessed, as well as the relationship between mathematical adjustment, simplicity and biological interpretation of parameters, according to the parsimony criterion (Santoyo et al. 2017).
Results and Discussion
Tables 2, 3 and 4 show the regression coefficients, as well as the components of variance of the error (σ2) and of the random effects (σ2 a). The goodness of fit criteria of models in the three mentioned approaches are exposed (without random effect, with random effect associated with the incubation bottle and with random effect linked to the bottle over time). The estimates obtained with the NLMIXED proc were those that maximized the integrated verisimilitude function on the random effects.
Logistic model |
Gompertz model |
|
||||||
---|---|---|---|---|---|---|---|---|
Fix parameter | ± SE | P Value | Fix parameter | ± SE | P Value | Fix parameter | ± SE | P Value |
b = 151.4 c = 0.04 L = 7.94 |
2.35 0.003 0.71 |
< 0.0001 < 0.0001 < 0.0001 |
b = 154.2 k = 5.27 r = 0.1 |
2 0.49 0.006 |
< 0.0001 < 0.0001 < 0.0001 |
b = 160 T = 20.4 F = 2.55 |
2.19 0.4 0.13 |
< 0.0001 < 0.0001 < 0.0001 |
σ2 | 61.98 | 39.6 | 26.1 | |||||
R2aj | 97.82 % | 98.61 % | 99.08 % | |||||
AIC | 285.6 | 268.2 | 249.9 | |||||
BIC | 297.3 | 279.8 | 259.9 | |||||
Residuals | ||||||||
DW | 1.0104 | 1.0161 | 0.8491 | |||||
K-S | p = 0.1378 | p = 0.0391 | p = 0.1500 | |||||
White | p = 0.6451 | p = 0.3457 | p = 0.0344 | |||||
Randomness | p = 0.009 | p = 0.023 | p = 0.023 |
The fermentation of the selected treatment (mixed silage + 25 % of sweet potato + L. pentosus LB-31) showed an asymptotic value of GP, which ranged between 138.9 and 160 mLg-1 incOM. GP rate was lower than 0.052 h-1. The constant factor of microbial efficiency was close to 0.1 (h-1). Lag phase was significant, indicating that GP did not start instantaneously. Half of GP was reached around 21 h. Similar results were obtained by Rodríguez et al. (2020) in a study of the same silage, but with the difference that the asymptotic GP values, estimated in this research, were close to 460 mL g-1 incOM.
The addition of random effects to the parameters favored the improvement of fit criteria, such as σ2, BIC and R2aj. Similar deductions were obtained by Gómez and Agüero (2020), when studying bovine lactation curves. Also, Corral et al. (2019) reached similar conclusions, when evaluating height/diameter relationship of seven Pinus species.
The R2aj were superior to 97 %, which indicates that models explained a large part of the observed variability. According to Rodríguez et al. (2019), sigmoid curves fit favorably to GP performance, since they allow describing the three phases of GP (rapid GP, deceleration, and linear or asymptotic phase). Out of the analyzed models, that of Groot et al. (1996) reached the most suitable values of R2aj, AIC and BIC. These results agree with those obtained by Chaves et al. (2021). These authors recommended the model of Groot et al. (1996) to estimate mean in vitro GP curves of silages. However, this model has been little known because it has a complex biological interpretation of the F parameter.
Regarding the basic assumptions for residuals, it was found that the residual mean was approximately zero in all cases. Normality and homogeneity of variances did not have a uniform performance in the evaluated models. Residual independence was not fulfilled. DW statistic was lower than 1.7, out of the required interval for the number of observations and parameters of the used models (Guerra et al. 2018). Autocorrelation, especially in the case of longitudinal data, could be due to the fact that variables that really affect the model have not been included or that the appropriate model has not been chosen (Gómez and Agüero 2020). Failure to comply with the assumption of independence of errors is a limitation that can lead to biased estimates (Pérez 2018).
The run test for randomness found that residuals were not random. The same was evident when studying residual graphs, in which a characteristic pattern was observed near 48 h. Once this happens, authors like Oddi et al. (2020) propose the inclusion of a variance function that allows modeling heteroscedasticity. However, no published information was found about specific functions of variance for in vitro GP curves.
In the second analysis, the incubation bottle effect was considered as a source of error that was not controlled in simple NLM. Table 3 presents the results. It was observed that the individual effects of bottles on the asymptotic GP were not statistically significant. Therefore, it was deduced that it was not necessary to add this effect in the regression analysis, and that, during the experiment, the incubation conditions of bottles were adequately controlled.
Logistic model |
Gompertz model |
Groot model |
||||||
---|---|---|---|---|---|---|---|---|
Fix parameter | ± SE | P Value | Fix parameter | ± SE | P Value | Fix parameter | ± SE | P Value |
b = 151.5 c = 0.04 L = 7.94 |
2.84 0.003 0.694 |
0.0004 0.004 0.008 |
b =154.2 k = 5.27 r = 0.1 |
2.8 0.454 0.005 |
0.0003 0.0074 0.0027 |
b =160 ta = 20.4 F = 2.54 |
3 0.35 0.11 |
0.0004 0.0003 0.0019 |
Random parameter | ± SE | P Value | Random parameter | ± SE | P Value | Random parameter | ± SE | P Value |
b1 = -2.31 b2 = 2.39 b3 = -0.09 |
3.46 3.53 2.43 |
0.5741 0.5676 0.9728 |
b1 = -3.48 b2 = 3.65 b3 = -0.17 |
3.01 3.05 2.64 |
0.37 0.35 0.95 |
b1 = -4.3 b2 = 4.59 b3 = -0.25 |
2.88 2.89 2.75 |
0.27 0.25 0.93 |
σ2 | 58.87 | 34.24 | 19.6 | |||||
σ2 a | 8.43 | 12.6 | 16.5 | |||||
R2aj | 98.02 % | 98.82 % | 99.33 % | |||||
AIC | 285.1 | 265.8 | 243.7 | |||||
BIC | 278.8 | 259.5 | 238.3 | |||||
Residuals | ||||||||
DW | 1.0897 | 1.1916 | 1.1143 | |||||
K-S | p = 0.15 | p = 0.0190 | p = 0.1500 | |||||
White | p = 0.5785 | p = 0.2533 | p = 0.0083 | |||||
Randomness | p = 0.009 | p = 0.105 | p = 0.023 |
Table 4 shows the results obtained, when the bottle-time random effect was added to explain the correlation between the repeated measures over time on the same experimental unit. It was observed that the growth of σ2 a caused a decrease in σ2. The induced correlation structure, when the bottle-time random effect was added, explained part of the variability of the GP data. With these results, it was evidenced that there was heterogeneity within bottles, and that it is important to consider this source of variation in the model. The bottle-time mixed NLMs reached the best fit, but did not solve the breach of the residual independence and randomness assumption. Although authors such as Corral et al. (2019) state that it is possible to minimize autocorrelation by adding random effects, the expected results were not obtained. The persistence of residual autocorrelation was attributed to the fact that the logistic, Gompertz and Groot et al. (1996) models have limitations in the description of in vitro GP.
Logistic model |
Gompertz model |
Groot model Groot |
||||||
---|---|---|---|---|---|---|---|---|
Fix parameter | ± SE | P Value | Fix parameter | ± SE | P Value | Fix parameter | ± SE | P Value |
b = 138.9 c = 0.052 L = 9.33 |
4.91 0.002 0.231 |
< 0.0001 < 0.0001 < 0.0001 |
b =154.1 k = 5.713 r = 0.106 |
2.49 0.242 0.003 |
0.0001 0.0001 0.0001 |
b =158.9 Ta = 20.2 F = 2.68 |
2.49 0.37 0.09 |
0.0001 0.0001 0.0001 |
σ2 | 1E-12 | 2.85 | 5.68 | |||||
σ2 a | 282.5 | 68.57 | 45.5 | |||||
R2aj | 99.99 % | 99.96 % | 99.91 % | |||||
AIC | 259.7 | 237.2 | 234.1 | |||||
BIC | 271.3 | 248.8 | 244.0 | |||||
Residuals | ||||||||
DW | 0.3183 | 1.2272 | 0.8468 | |||||
K-S | p = 0.01 | p = 0.0100 | p = 0.0159 | |||||
White | p = 0.0003 | p = 0.1424 | p = 0.1109 | |||||
Randomness | p = 0.001 | p = 0.004 | p = 0.023 |
The estimates of the random effects of each model can be observed in tables 5, 6 and 7. It is important to note that each set of bij behaved in a normal way, with zero mean and constant variance. The addition of the bij effect to the models produced a characteristic GP for each bottle, according to the time. In addition, differences between the asymptotic GP b + bij and the general mean b were reflected (Bandera and Pérez 2018). Through this analysis, it was observed that the 48 h of incubation was decisive in the estimation of the in vitro GP asymptote. At approximately this time, the NLM partially characterized the GP. This could be due to limitations of the models used to adequately describe in vitro GP in all its stages or because the readings after the first 24 h were not equally spaced. By not measuring GP frequently and systematically, fluctuations typical to the in vitro GP process may no longer be recorded (Muro et al. 2017).
Time | i =1 | ± SE | P value | i =2 | ± SE | P value | i = 3 | ± SE | P value |
---|---|---|---|---|---|---|---|---|---|
j = 3 h | -14.72 | 6.46 | 0.028 | 3.90 | 7.34 | 0.598 | 15.94 | 7.94 | 0.052 |
j = 6 h | -14.78 | 4.86 | 0.004 | 1.01 | 5.36 | 0.851 | -23.28 | 4.63 | 0.000 |
j = 9 h | -8.10 | 4.38 | 0.072 | 6.04 | 4.72 | 0.209 | -7.41 | 4.39 | 0.100 |
j = 12 h | 14.59 | 5.39 | 0.010 | 28.43 | 5.81 | 0.000 | 16.22 | 5.44 | 0.005 |
j = 15 h | 15.16 | 6.26 | 0.020 | 26.25 | 6.62 | 0.000 | 21.79 | 6.48 | 0.002 |
j = 18 h | 5.73 | 6.38 | 0.375 | 12.47 | 6.59 | 0.066 | 14.03 | 6.64 | 0.041 |
j = 21 h | -0.99 | 6.03 | 0.871 | 7.84 | 6.25 | 0.218 | 8.99 | 6.28 | 0.160 |
j = 24 h | -7.43 | 5.48 | 0.183 | -0.44 | 5.59 | 0.938 | 2.40 | 5.63 | 0.673 |
j = 29h | -16.73 | 4.93 | 0.002 | -7.53 | 4.97 | 0.138 | -6.78 | 4.97 | 0.180 |
j = 48h | -10.49 | 4.90 | 0.039 | -0.02 | 4.90 | 0.997 | -3.04 | 4.90 | 0.538 |
j = 72h | 7.68 | 4.91 | 0.126 | 18.62 | 4.91 | 0.001 | 10.69 | 4.91 | 0.036 |
j = 77h | 14.13 | 4.91 | 0.007 | 25.19 | 4.91 | 0.000 | 15.10 | 4.91 | 0.004 |
j = 144h | 22.14 | 4.91 | 0.000 | 33.35 | 4.91 | 0.000 | 21.48 | 4.91 | 0.000 |
Time | i=1 | ± SE | P value | i=2 | ± SE | P value | i=3 | ± SE | P value |
---|---|---|---|---|---|---|---|---|---|
j = 3h | 0.74 | 8.27 | 0.929 | 0.99 | 8.27 | 0.906 | 1.15 | 8.28 | 0.891 |
j = 6h | 0.48 | 8.08 | 0.953 | 1.59 | 8.11 | 0.846 | -0.12 | 8.08 | 0.988 |
j = 9h | -4.86 | 7.76 | 0.535 | -1.6 | 7.47 | 0.832 | -4.7 | 7.74 | 0.548 |
j = 12h | -4.54 | 6.27 | 0.473 | 1.91 | 6.18 | 0.759 | -3.78 | 6.23 | 0.548 |
j = 15h | -2.76 | 4.89 | 0.576 | 4.8 | 5.04 | 0.346 | 1.76 | 4.94 | 0.723 |
j = 18h | -2.32 | 4.09 | 0.573 | 3.42 | 4.2 | 0.421 | 4.75 | 4.24 | 0.270 |
j = 21h | -0.36 | 3.64 | 0.921 | 8.23 | 3.83 | 0.038 | 9.36 | 3.87 | 0.020 |
j = 24h | -2.05 | 3.26 | 0.534 | 5.27 | 3.38 | 0.127 | 8.24 | 3.44 | 0.022 |
j = 29h | -11.97 | 2.81 | 0.000 | -2.04 | 2.85 | 0.478 | -1.24 | 2.85 | 0.667 |
j = 48h | -20.39 | 2.77 | 0.000 | -10.03 | 2.7 | 0.001 | -13.03 | 2.72 | 0.000 |
j = 72h | -6.74 | 2.89 | 0.025 | 3.79 | 2.86 | 0.193 | -3.84 | 2.88 | 0.191 |
j = 77h | -0.69 | 2.89 | 0.812 | 9.94 | 2.86 | 0.001 | 0.24 | 2.88 | 0.935 |
j = 144h | 6.75 | 2.89 | 0.025 | 17.52 | 2.89 | 0.000 | 6.12 | 2.9 | 0.041 |
Time | i=1 | ± SE | P value | i=2 | ± SE | P value | I = 3 | ± SE | P value |
---|---|---|---|---|---|---|---|---|---|
j = 3h | 0.16 | 6.74 | 0.981 | 0.2 | 6.74 | 0.977 | 0.22 | 6.74 | 0.975 |
j = 6h | 0.59 | 6.71 | 0.930 | 0.89 | 6.72 | 0.896 | 0.43 | 6.71 | 0.949 |
j = 9h | -1.24 | 6.54 | 0.851 | -0.04 | 6.51 | 0.995 | -1.18 | 6.54 | 0.858 |
j = 12h | -2.82 | 6.1 | 0.647 | 0.36 | 5.99 | 0.953 | -2.44 | 6.07 | 0.690 |
j = 15h | -3.55 | 5.31 | 0.508 | 1.17 | 5.25 | 0.825 | -0.73 | 5.23 | 0.889 |
j = 18h | -3.5 | 4.58 | 0.449 | 0.7 | 4.56 | 0.879 | 1.67 | 4.58 | 0.718 |
j = 21h | -0.91 | 4.08 | 0.826 | 6.03 | 4.28 | 0.167 | 6.94 | 4.33 | 0.117 |
j = 24h | -0.73 | 3.78 | 0.848 | 5.57 | 3.91 | 0.163 | 8.12 | 4.02 | 0.050 |
j = 29h | -7.52 | 3.45 | 0.036 | 1.57 | 3.44 | 0.650 | 2.31 | 3.45 | 0.508 |
j = 48h | -15.72 | 3.13 | 0.000 | -5.76 | 2.93 | 0.057 | -8.64 | 2.97 | 0.006 |
j = 72h | -6.48 | 3 | 0.037 | 3.5 | 2.93 | 0.241 | -3.73 | 2.97 | 0.216 |
j = 77h | -1.33 | 2.97 | 0.656 | 8.71 | 2.97 | 0.006 | -0.46 | 2.96 | 0.879 |
j = 144h | 2.74 | 3.09 | 0.380 | 12.75 | 3.14 | 0.000 | 2.16 | 3.09 | 0.489 |
It is concluded that the non-linear mixed models allowed to obtain better R2aj, BIC and σ2. Random effects did not have a significant effect on the fulfillment of the assumptions of autocorrelation and randomness of residuals. The random effect bottle over time allowed to determine the times in which there was a lack of fit. Despite limitations of the analyzed models, the logistic, without the inclusion of the random effect, was the most appropriate according to the parsimony criterion.