Acessibilidade / Reportar erro

Bayesian modeling of the coffee tree growth curve

Modelagem bayesiana da curva de crescimento do cafeeiro

ABSTRACT:

When modeling growth curves, it should be considered that longitudinal data may show residual autocorrelation, and, if this characteristic is not considered, the results and inferences may be compromised. The Bayesian approach, which considers priori information about studied phenomenon has been shown to be efficient in estimating parameters. However, as it is generally not possible to obtain marginal distributions analytically, it is necessary to use some method, such as the weighted resampling method, to generate samples of these distributions and thus obtain an approximation. Among the advantages of this method, stand out the generation of independent samples and the fact that it is not necessary to evaluate convergence. In this context, the objective of this work research was: to present the Bayesian nonlinear modeling of the coffee tree height growth, irrigated and non-irrigated (NI), considering the residual autocorrelation and the nonlinear Logistic, Brody, von Bertalanffy and Richard models. Among the results, it was found that, for NI plants, the Deviance Information Criterion (DIC) and the Criterion of density Predictive Ordered (CPO), indicated that, among the evaluated models, the Logistic model is the one that best describes the height growth of the coffee tree over time. For irrigated plants, these same criteria indicated the Brody model. Thus, the growth of the non-irrigated and irrigated coffee tree followed different growth patterns, the height of the non-irrigated coffee tree showed sigmoidal growth with maximum growth rate at 726 days after planting and the irrigated coffee tree starts its development with high growth rates that gradually decrease over time.

Key words:
residual autocorrelation; nonlinear models; Logistic model; Brody model; Von Bertalanffy model; Richards model.

RESUMO:

Na modelagem de curvas de crescimento deve-se considerar que dados longitudinais podem apresentar autocorrelação residual, sendo que, se tal característica não é considerada, os resultados e inferências podem ser comprometidos. A abordagem bayesiana, que considera informações à priori sobre o fenômeno em estudo tem se mostrado eficiente na estimação de parâmetros. No entanto, como geralmente não é possível obter as distribuições marginais de forma analítica, faz-se necessário a utilização de algum método, como o método de reamostragem ponderada, para gerar amostras dessas distribuições e assim obter uma aproximação para as mesmas. Dentre as vantagens desse método, destaca-se a geração de amostras independentes e o fato de não ser necessário avaliar convergência. Diante desse contexto, o objetivo deste trabalho foi apresentar a modelagem não linear bayesiana do crescimento em altura de plantas do cafeeiro, irrigadas e não irrigadas (NI), considerando a autocorrelação residual e os modelos não lineares Logístico, Brody, von Bertalanffy e Richards. Em vista dos resultados, verificou-se que, para as plantas NI, o DIC e CPOc, indicaram que, dentre os modelos avaliados, o modelo Logístico é o que melhor descreve o crescimento em altura do cafeeiro ao longo do tempo. E, para as plantas irrigadas, esses mesmos critérios indicaram o modelo Brody. Assim, o crescimento da planta do cafeeiro não irrigado e irrigado seguiram padrões de crescimento distintos, a altura do cafeeiro não irrigado apresentou crescimento sigmoidal com taxa máxima de crescimento aos 726 dias após o plantio, já o cafeeiro irrigado inicia seu desenvolvimento com altas taxas de crescimento que vão diminuindo aos poucos com o tempo.

Palavras-chave:
autocorrelação residual; modelos não lineares; modelo Logístico; modelo Brody; modelo Von Bertalanffy; Modelo Richards

INTRODUCTION:

Due to the economic importance of coffee for Brazil, one of the main focuses of research in the agricultural sector refers to the development of coffee tree. A relevant variable to be considered in such research is the plant height, because in addition to being representative of vegetative development, it is also correlated with productivity (ASSIS et al., 2014ASSIS, G. A. de et al. Correlation between coffee plant growth and yield as function of water supply regime and planting density. . Bioscience Journal, Uberlândia, v.30, n.3, p.666-676, 2014. Available from: <Available from: http://www.seer.ufu.br/index.php/biosciencejournal/article/view/18044 >. Accessed: Feb. 01, 2021.
http://www.seer.ufu.br/index.php/bioscie...
; PEREIRA et al., 2016PEREIRA, A. A. et al. Modeling nonlinear growth in height coffee with and without irrigation in different densities. Irriga, Botucatu, v.1, n.1, p.140-149, 2016. Available from: <Available from: http://irriga.fca.unesp.br/index.php/irriga/article/view/1863 >. Accessed: Feb. 01, 2021. doi: https://doi.org/10.15809/irriga.2016v1n1p140-149.
http://irriga.fca.unesp.br/index.php/irr...
; SILVA et al., 2019SILVA, E. M. et al. Leaf count overdispersion in coffee seedlings. Ciência Rural, Santa Maria, v.49, n.4, p.1-7, 2019. Available from: <Available from: https://www.scielo.br/scielo.php?script=sci_arttext&pid=S0103-84782019000400201 >. Accessed: Feb. 01, 2021. doi: 10.1590/0103-8478cr20180786.
https://www.scielo.br/scielo.php?script=...
).

Modeling of the coffee tree growth curve using nonlinear regression models allows the researcher to describe growth over time and obtain information of practical interest, such as the maximum height at adult stage and the speed of growth. They are important to characterize the development of plants (PEREIRA et al., 2014PEREIRA, A. A. et al. Description vegetative growth of coffee tree farming ruby mg 1192 using regression models. Coffee Science, Lavras, v.9, n.2, p.266-274, abr./jun. 2014. Available from: <http://www.sDICafe.ufv.br/handle/123456789/8051>. Accessed: Feb. 01, 2021. ; REIS et al., 2014REIS, R. M. et al. Nonlinear regression models applied to clusters of garlic accessions. Horticultura Brasileira, Vitoria da Conquista, v.32, n.2, p.178-183, abr./jun., 2014. Available from: <Available from: https://www.scielo.br/scielo.php?pid=S0102-05362014000200178&script=sci_abstract&tlng=pt >. Accessed: Feb. 01, 2021. doi: 10.1590/S0102-05362014000200010.
https://www.scielo.br/scielo.php?pid=S01...
; FERNANDES et al., 2017FERNANDES, T. J. et al. Double sigmoidal models describing the growth of coffee berries. Ciência Rural, Santa Maria, v.47, n.8, p.1-7, 2017. Available from: <Available from: https://www.scielo.br/scielo.php?script=sci_arttext&pid=S0103-84782017000800401 >. Accessed: Feb. 01, 2021. doi: 10.1590/0103-8478cr20160646.
https://www.scielo.br/scielo.php?script=...
; SARI et al., 2018SARI, B. G. et al. Nonlinear modeling for analyzing data from multiple harvest crops. Agronomy Journal, v.110, n.6, p.2331-2342, 2018. Available from: <Available from: https://acsess.onlinelibrary.wiley.com/doi/full/10.2134/agronj2018.05.0307 >. Accessed: Feb. 01, 2021. doi: 10.2134/agronj2018.05.0307
https://acsess.onlinelibrary.wiley.com/d...
; SILVA & SAVIAN, 2019SILVA, P. V. & SAVIAN, T. V. Chanter model: nonlinear modeling of the fruit growth of cocoa. Ciência Rural, Santa Maria, v.49, n.11, p.1-7, 2019. Available from: <Available from: https://www.scielo.br/scielo.php?script=sci_arttext&pid=S0103-84782019001100407 >. Accessed: Feb. 01, 2021. doi: 10.1590/0103-8478cr20190409
https://www.scielo.br/scielo.php?script=...
; JANE et al., 2020JANE, S. A. et al. Adjusting the growth curve of sugarcane varieties using nonlinear models. Ciência Rural, Santa Maria, v.50, n.3, p.1-7, 2020. Available from: <Available from: https://www.scielo.br/scielo.php?pid=S0103-84782020000300204&script=sci_arttext >. Accessed: Feb. 01, 2021. doi: 10.1590/0103-8478cr20190408.
https://www.scielo.br/scielo.php?pid=S01...
; SILVA et al., 2021; MENDES et al., 2021MENDES, P. N. et al. Logistic Bayesian model in the study of the growth of tomatoes . Research, Society and Development, v.10, n.3, p.1-12, 2021. Available from: <Available from: https://rsdjournal.org/index.php/rsd/article/view/13198 >. Accessed: Apr. 05, 2021. doi: 10.33448/rsd-v10i3.13198
https://rsdjournal.org/index.php/rsd/art...
).

To model the growth curve, it should be considered that the growth data taken over time may show residual autocorrelation (PEREIRA et al., 2016PEREIRA, A. A. et al. Modeling nonlinear growth in height coffee with and without irrigation in different densities. Irriga, Botucatu, v.1, n.1, p.140-149, 2016. Available from: <Available from: http://irriga.fca.unesp.br/index.php/irriga/article/view/1863 >. Accessed: Feb. 01, 2021. doi: https://doi.org/10.15809/irriga.2016v1n1p140-149.
http://irriga.fca.unesp.br/index.php/irr...
; JANE et al., 2020JANE, S. A. et al. Adjusting the growth curve of sugarcane varieties using nonlinear models. Ciência Rural, Santa Maria, v.50, n.3, p.1-7, 2020. Available from: <Available from: https://www.scielo.br/scielo.php?pid=S0103-84782020000300204&script=sci_arttext >. Accessed: Feb. 01, 2021. doi: 10.1590/0103-8478cr20190408.
https://www.scielo.br/scielo.php?pid=S01...
; SILVA et al., 2021SILVA, W. S. et al. Eucalyptus grandis x eucalyptus urophylla growth curve in different site classifications, considering residual autocorrelation. Revista Brasileira de Biometria, v.39, n.1, p.122-138, 2021. Available from: <Available from: https://biometria.ufla.br/index.php/BBJ/article/view/511 >. Accessed: Jun. 01, 2021. doi: 10.28951/rbb.v39i1.511
https://biometria.ufla.br/index.php/BBJ/...
). However, in several applied researches, it appears that this characteristic is not incorporated in the modeling, and may compromise the results and inferences (MUNIZ et al., 2017MUNIZ, J. A. et al. Nonlinear models for description of cacao fruit growth with assumption violations. Caatinga, Mossoró, v.30, n.1, p.250-257, jan./mar. 2017. Available from: <Available from: https://www.scielo.br/scielo.php?script=sci_arttext&pid=S1983-21252017000100250 >. Accessed: Feb. 01, 2021. doi: 10.1590/1983-21252017v30n128rc.
https://www.scielo.br/scielo.php?script=...
).

An approach that has been shown to be efficient in estimating parameters is the Bayesian approach, which considers both observations and model parameters as random variables. Thus, an a prior distribution for the parameters is specified, denoted by P (θ) - which represents in terms of probability the pre-existing knowledge about the parameters - and a joint probability distribution for the sampling data, called the likelihood function L (θ ǀ Y) (BOX & TIAO, 1992BOX, G. E. P. & TIAO, G. C. Bayesian inference in statistical analysis. New York: J. Wiley,1992. 588p. ; PAULINO et al., 2018PAULINO, C. D. et al. Estatística Bayesiana. Lisboa: Fundação Calouste Gulbenkian, 2º edição, 2018. 601p.; BOLSTAD & CURRAN, 2016BOLSTAD, W. M. & CURRAN, J. M. Introduction to Bayesian Statistics. Third Edition, 2016.). The combination of this information, by Bayes’ theorem, results in the posterior distribution of θ. This distribution, indicated by P (θ ǀ Y) represents the updated knowledge about the parameters, associating the researcher initial knowledge with the information from the sample: Pθ| Y=Lθ| YPθPY=Lθ| YPθPθ, YdθLθ| YPθ, where P(Y) is a normalizing constant.

Among the advantages of the Bayesian approach, stands out satisfactory modeling even with a relatively small sample and obtaining credibility intervals (MARTINS FILHO et al., 2008MARTINS FILHO, S. et al. Bayesian approach in the growth curves of two cultivars of common bean. Ciência Rural, Santa Maria, v.38, n.6, p.1516-1521, 2008. Available from: <Available from: https://www.scielo.br/scielo.php?pid=S0103-84782008000600004&script=sci_abstract&tlng=pt >. Accessed: Feb. 01, 2021. doi: 10.1590/S0103-84782008000600004.
https://www.scielo.br/scielo.php?pid=S01...
; SILVA et al., 2020SILVA, E. M. et al. Bayesian approach to the zinc extraction curve of soil with sewage sludge. Acta Scientiarum. Technology, v.42, n.1, p.1-9, 2020. Available from: <Available from: http://periodicos.uem.br/ojs/index.php/ActaSciTechnol/article/view/46893 >. Accessed: Feb. 01, 2021. 10.4025/actascitechnol.v42i1.46893
http://periodicos.uem.br/ojs/index.php/A...
; MACEDO et al., 2017MACEDO, L. R. et al. Bayesian inference for the fitting of dry matter accumulation curves in garlic plants. Pesquisa Agropecuária Brasileira, v.52, n.8, p.572-581, 2017. Available from: <Available from: https://www.scielo.br/scielo.php?script=sci_arttext&pid=S0100-204X2017000800572 >. Accessed: Feb. 01, 2021. doi: 10.1590/s0100-204x2017000800002.
https://www.scielo.br/scielo.php?script=...
).

In order to make inferences about a given parameter θi, i= 1, 2, ..., p, the posterior marginal distribution of this parameter should be obtained, which is given by the integral of the posterior distribution in relation to the other parameters. Due to the complexity of these integrals, in most cases, obtaining the marginal distributions is not done in an analytical way, making it necessary to use some approximation method. Markov Chain Monte Carlo methods (MCMC) have been the most used in the generation of samples from a posterior distribution, thus making it possible to obtain an approximate density. However, as they are iterative, these methods may take time to converge and require the researcher to use convergence assessment diagnostics (COLE et al., 2012COLE, S. R. et al. Bayesian posterior distributions without markov chains. American Journal of Epidemiology, Oxford, v.175, n.5, p.368-375, Feb. 2012. Available from: <Available from: https://pubmed.ncbi.nlm.nih.gov/22306565/ >. Accessed: Feb. 1, 2021.
https://pubmed.ncbi.nlm.nih.gov/22306565...
).

An alternative method to the MCMC, which generates independent samples, does not require convergence assessment and is computationally faster, is the weighted resampling method, proposed by RUBIN (1987RUBIN, D. B. A noniterative sampling/importance resampling alternative to the data augmentation algorithm for creating a few imputations when fractions of missing information are modest: the SIR algorithm. Journal of the American Statistical Association, New York, v.82, n.398, p.543-546, 1987. ) (COLE et al., 2012COLE, S. R. et al. Bayesian posterior distributions without markov chains. American Journal of Epidemiology, Oxford, v.175, n.5, p.368-375, Feb. 2012. Available from: <Available from: https://pubmed.ncbi.nlm.nih.gov/22306565/ >. Accessed: Feb. 1, 2021.
https://pubmed.ncbi.nlm.nih.gov/22306565...
; LOPES et al., 2012LOPES, H. F. et al. Bayesian statistics with a smile: a resampling-sampling perspective. Brazilian Journal of Probability and Statistics, São Paulo, v.26, n.4, p.358-371, 2012. Available from: <Available from: https://www.jstor.org/stable/43601224?seq=1 >. Accessed: Feb. 01, 2021.
https://www.jstor.org/stable/43601224?se...
). The generation of samples by this method, also known as sampling-importance resampling (SIR), is done in two stages: initially, the values of a candidate distribution are generated, and in the second stage, resampling is performed considering the weight (probability) assigned to each value generated. According with COLE et al. (2012COLE, S. R. et al. Bayesian posterior distributions without markov chains. American Journal of Epidemiology, Oxford, v.175, n.5, p.368-375, Feb. 2012. Available from: <Available from: https://pubmed.ncbi.nlm.nih.gov/22306565/ >. Accessed: Feb. 1, 2021.
https://pubmed.ncbi.nlm.nih.gov/22306565...
) the SIR has a more attractive and easy-to-implement theory, so that researchers can focus their efforts on eliciting prior distributions.

Although, there are applied researches in the literature presenting the fitting of regression models using Bayesian methods, such researches generally use linear and nonlinear models . They are limited to the assumption of residual independence (ANDRADE FILHO et al.; 2010ANDRADE FILHO, M. G. de et al. Bayesian approach to growth curve with objective priori. Revista Brasileira de Biometria, São Paulo, v.28, n.2, p.161-181, 2010. Available from: <Available from: http://jaguar.fcav.unesp.br/RME/fasciculos/v28/v28_n2/A10_Marinho.pdf >. Accessed: Feb. 01, 2021.
http://jaguar.fcav.unesp.br/RME/fascicul...
; MARTINS FILHO et al.; 2008MARTINS FILHO, S. et al. Bayesian approach in the growth curves of two cultivars of common bean. Ciência Rural, Santa Maria, v.38, n.6, p.1516-1521, 2008. Available from: <Available from: https://www.scielo.br/scielo.php?pid=S0103-84782008000600004&script=sci_abstract&tlng=pt >. Accessed: Feb. 01, 2021. doi: 10.1590/S0103-84782008000600004.
https://www.scielo.br/scielo.php?pid=S01...
; SILVA et al., 2020SILVA, E. M. et al. Bayesian approach to the zinc extraction curve of soil with sewage sludge. Acta Scientiarum. Technology, v.42, n.1, p.1-9, 2020. Available from: <Available from: http://periodicos.uem.br/ojs/index.php/ActaSciTechnol/article/view/46893 >. Accessed: Feb. 01, 2021. 10.4025/actascitechnol.v42i1.46893
http://periodicos.uem.br/ojs/index.php/A...
; MACEDO et al., 2017MACEDO, L. R. et al. Bayesian inference for the fitting of dry matter accumulation curves in garlic plants. Pesquisa Agropecuária Brasileira, v.52, n.8, p.572-581, 2017. Available from: <Available from: https://www.scielo.br/scielo.php?script=sci_arttext&pid=S0100-204X2017000800572 >. Accessed: Feb. 01, 2021. doi: 10.1590/s0100-204x2017000800002.
https://www.scielo.br/scielo.php?script=...
) or consider autocorrelation and use linear models (CHIB & GREENBERG, 1994CHIB, S. & GREENBERG, E. Bayes inference in regression models with ARMA (p, q) errors. Journal of Econometrics, North Holland, v.64, p.183-206, 1994. Available from: <Available from: https://www.sciencedirect.com/science/article/abs/pii/0304407694900639 >. Accessed: Feb. 01, 2021. doi: 10.1016/0304-4076(94)90063-9.
https://www.sciencedirect.com/science/ar...
; MENZEFRICKE, 1999MENZEFRICKE, U. Bayesian prediction in growth-curve models with correlated errors. Test, Pamplona, v.8, n.1, p.75-93, 1999. Available from: <Available from: https://link.springer.com/article/10.1007/BF02595863 >. Accessed: Feb. 01, 2021.
https://link.springer.com/article/10.100...
). Specifically, in the study of coffee tree growth, there is no research in the literature using a Bayesian approach considering both nonlinear models and residual autocorrelation.

Given this context, the objective of this study was: to present the Bayesian nonlinear modeling of coffee tree height growth, irrigated and non-irrigated, considering residual autocorrelation and nonlinear Logistic, Brody, von Bertalanffy and Richards models.

MATERIALS AND METHODS:

Data analyzed comes from an experiment carried out with the cultivar Rubi MG1192, in the experimental area of the Department of Agriculture, Federal University of Lavras, in Lavras, state of Minas Gerais.

Coffee crop was planted in January 2001, using the 4m x 1m spacing, which corresponds to the planting of 2,500 plants per hectare. During the conduct of the experiment, crop treatments and phytosanitary control were carried out in accordance with the requirements of the crop. Fertilizer application and liming were done according to the analysis of soil and micronutrients provided by foliar fertilization.

The design used was in randomized blocks, with four replications. In each block, six irrigation regimes were drawn; however, in this study, only two regimes were analyzed: non-irrigated (NI) and irrigated at 60 kPa (most commonly used in this region). The irrigation system implanted in the field was the drip system. Irrigations at 60 kPa occurred when the tensiometer, at a depth of 25 cm, recorded this tension value on the tensimeter, which is measured in kilopascal (kPa).

Each experimental unit consisted of three planting rows with ten plants in each row, the eight central plants of the central row were considered useful. The plant height was represented by the average of the eight useful plants measured in centimeters (cm), measured quarterly, in the period between May 2001 and August 2006, totaling 22 measurements.

To analyze this growth, the nonlinear regression models in table 1 were used (PEREIRA et al., 2014PEREIRA, A. A. et al. Description vegetative growth of coffee tree farming ruby mg 1192 using regression models. Coffee Science, Lavras, v.9, n.2, p.266-274, abr./jun. 2014. Available from: <http://www.sDICafe.ufv.br/handle/123456789/8051>. Accessed: Feb. 01, 2021. ; PEREIRA et al., 2016). Considering that the residuals (white noise) follow a Normal distribution, that is, ui ~ N(0, σ2 u), we have the following expression for the likelihood function: LθY=i=12212πσu2exp-1(yi-μi)22σu2..

Table 1
Expected value and autocorrelation structure for the residuals of the Logistic, Brody, von Bertalanffy and Richards models considering the data referring to the non-irrigated (NI) and irrigated regimes.

Y represents the vector formed by the observed values yi and, the expected value μi and θ vary according to the model and the presence and order of residual autocorrelation, which were verified through the graphs of the autocorrelation (acf) and partial autocorrelation (pacf) functions of the residuals of the ordinary model (considering all assumptions met). In cases where acf suggested residual dependency, pacf was used to identify the p order of the AR (if AR(1), AR(2), etc) (MORETTIN & TOLOI, 2006MORETTIN, P. A. & TOLOI, C. M. C. Análise de séries temporais. 2. ed. São Paulo: Blucher, 2006. 538 p.).

It should be noted that the equations presented in table 1 represent the expression for the mean µ that will be substituted in the likelihood presented in the previous paragraph. In this parameterizations, ti indicated the evaluation times, in days after planting: ti = 90, 180, ..., 1980 days; α is the upper horizontal asymptote, that is, it represents the maximum expected height; β indicates the abscissa of the inflection point in the Logistical and Richards models, in Brody and Von Bertalanffy it is a scale parameter without direct biological interpretation; δ is the parameter that defines the shape of the curve and determines in what proportion of α the inflection point occurs; k indicates the growth rate and determines the efficiency of growth: the higher this value, the less time it will take for α to be reached; ei is the residual which is assumed to be independently and identically distributed according to a Normal distribution with zero mean and homogeneous variance, that is, ei ~ N(0, σ2).

Thus, the adjustments that presented residual autocorrelation of order 1 were redone by the generalized least squares method, incorporating in the modeling the autoregressive of order 1, φ1 parameter, given by: ei= φ1ei-1+ui.

Similarly, for those who presented autocorrelation of order 2, the autoregressive parameter of order 2, φ2 was incorporated:

ei=φ1ei-1+φ2ei-2+ui,

Where ei, ei-1 e ei-2 correspond to the residuals generated in the times t i , t i - 1 e t i - 2 respectively and uiis the white noise, obviously in a structure without autocorrelation ei=ui (MORETTIN & TOLOI, 2006MORETTIN, P. A. & TOLOI, C. M. C. Análise de séries temporais. 2. ed. São Paulo: Blucher, 2006. 538 p.).

The prior distributions were elicited based on the technique known as “prior of specialist”, which determined a plausible range of occurrence for the parameters and some measure of position, such as the mean, according to an expert’s knowledge of the phenomenon studied and about parametric space of parameters. According to MOALA & PENHA (2016MOALA, F. A. & PENHA, D. L. Elicitation methods for beta prior distribution. Revista Brasileira de Biometria, Lavras, v.34, n.1, p.49-62, jan. 2016. Available from: <Available from: http://www.biometria.ufla.br/index.php/BBJ/article/view/91 >. Accessed: Feb. 01, 2021.
http://www.biometria.ufla.br/index.php/B...
), the elicitation of a prior distributions for the parameters constitutes an important phase in carrying out the Bayesian analysis. Thus, after exploratory analysis of the behavior of the growth curve parameters and about Gamma and Beta distributions, the following priorities were elicited.

For parameter α, which indicated an estimate of the maximum height to be reached by coffee trees; therefore α > 0 (SARI et al., 2018SARI, B. G. et al. Nonlinear modeling for analyzing data from multiple harvest crops. Agronomy Journal, v.110, n.6, p.2331-2342, 2018. Available from: <Available from: https://acsess.onlinelibrary.wiley.com/doi/full/10.2134/agronj2018.05.0307 >. Accessed: Feb. 01, 2021. doi: 10.2134/agronj2018.05.0307
https://acsess.onlinelibrary.wiley.com/d...
), a Gamma distribution with hyperparameters 10 and 0.04 was considered as a prior distribution.

For parameter β, which in the Logistic and Richards models, indicated the moment when the growth rate is maximum (inflection point) and in the other models does not have a direct practical interpretation (SARI et al., 2018SARI, B. G. et al. Nonlinear modeling for analyzing data from multiple harvest crops. Agronomy Journal, v.110, n.6, p.2331-2342, 2018. Available from: <Available from: https://acsess.onlinelibrary.wiley.com/doi/full/10.2134/agronj2018.05.0307 >. Accessed: Feb. 01, 2021. doi: 10.2134/agronj2018.05.0307
https://acsess.onlinelibrary.wiley.com/d...
), a Gamma distribution with hyperparameters 1.5 and 0.003 was also considered.

Parameter k, which is interpreted as a growth index, generally assumes low values, concentrated near the lower limit of the interval (0.1) (SARI et al., 2018SARI, B. G. et al. Nonlinear modeling for analyzing data from multiple harvest crops. Agronomy Journal, v.110, n.6, p.2331-2342, 2018. Available from: <Available from: https://acsess.onlinelibrary.wiley.com/doi/full/10.2134/agronj2018.05.0307 >. Accessed: Feb. 01, 2021. doi: 10.2134/agronj2018.05.0307
https://acsess.onlinelibrary.wiley.com/d...
). Thus, a Beta distribution (2.25; 1,500) was considered as a prior distribution.

When using the Richards model, specifying an a prior distribution for the δ parameter is also necessary. For this parameter, which can assume positive or negative values, a Normal was considered with a mean of 0.11 and standard deviation 10.

A generalized Beta distribution was defined as a prior for φ1 whose occurrence interval is [-1, 1], because according to MORETTIN & TOLOI (2006MORETTIN, P. A. & TOLOI, C. M. C. Análise de séries temporais. 2. ed. São Paulo: Blucher, 2006. 538 p.), for the process to be stationary, the following condition must be satisfied: |φ1|<1.

The generalized Beta distribution is obtained as follows: whether V ~ Beta (h1, h2) defined in the interval (0, 1), then the random variable defined by W= (n-m)V+m in the finite interval [m, n] follows a generalized Beta distribution (h1, h2) (MCDONALD & XU, 1995MCDONALD, J. B. & XU, Y. J. A generalization of the beta distribution with applications. Journal of Econometrics, Amsterdam, v.66, n.1, p.133-152, Mar./Apr. 1995. Available from: <Available from: https://www.sciencedirect.com/science/article/abs/pii/0304407694016124#:~:text=A%20generalization%20of%20the%20beta%20distribution%20with%20applications%E2%98%86&text=The%20generalized%20beta%20leads%20to,the%20normal%20as%20special%20cases >. Accessed: Feb. 01, 2021. doi: 10.1016/0304-4076(94)01612-4.
https://www.sciencedirect.com/science/ar...
). Thus, when considering V ~ Beta (3, 2), the random variable φ1 defined as φ1 = (1 + 1).V-1 has a generalized Beta distribution (3, 2) defined in the interval [-1, 1]. Therefore, the prior distribution of φ1 is given by [2.Beta (3, 2) -1].

The prior distribution established for parameter φ2 was a Uniform (-1, 1), with the following stationary conditions being observed: φ1+ φ2<1, φ2 - φ1<1 and | φ2|<1.

Hyperparameters determination of the prior distributions (which was done in order to obtain large variance) and the attribution of an non-informative distribution for parameter φ2 is based on the fact that the literature is scarce in research in the plant area related to Bayesian growth curve modeling considering the autocorrelation between the residuals.

For each model analyzed, independence between the parameters was considered. Thus, the joint prior distribution for the vector of parameters of interest θ, P(θ), is given by the prior product of each parameter. The number of parameters in θ varies depending on the model and the autocorrelation structure indicated by the graphic of pacf.

Using Bayes’ theorem and replacing the expressions obtained for the likelihood function and the joint a prior distribution, and discarding the integration constant, we have that the a posterior distribution P (θ | Y), for each case, is given by:

PθYexp-12σu2i=122yi-μi2P(θ).

To obtain an informative summary on the a posterior marginal distribution of each parameter, it is necessary to solve integrals of the form: EhθY=hθYP(θ|Y)dθ, where h(θ) is a conveniently chosen function. However, as it is not possible to solve these integrals analytically, the sampling-importance resampling (SIR) based on SMITH & GELFAND (1992SMITH A. F. M., & GELFAND A.E. Bayesian statistics without tears: a sampling-resampling perspective. The American Statistician. v.46, p.84-88, 1992. Available from: <Available from: https://www.tandfonline.com/doi/abs/10.1080/00031305.1992.10475856 >. Accessed Jan. 23, 2021. doi: 10.1080/00031305.1992.10475856
https://www.tandfonline.com/doi/abs/10.1...
) theory, was used to generate samples. Subsequently, estimates of the mean, mode and highest a posterior density (HPD) were obtained for the parameters (PAULINO et al., 2018PAULINO, C. D. et al. Estatística Bayesiana. Lisboa: Fundação Calouste Gulbenkian, 2º edição, 2018. 601p.).

To indicate the model that presented the best fit, considering the Bayesian approach, the Deviance information criterion - DIC (SPIEGELHALTER et al., 2002SPIEGELHALTER, D. J. et al. Bayesian measures of model complexity and fit. Journal of the Royal Statistical Society: Series B, v.64, n.4, p.583-639, 2002. Available from: <Available from: https://rss.onlinelibrary.wiley.com/doi/10.1111/1467-9868.00353 >. Accessed: Jun. 01, 2021. doi: 10.1111/1467-9868.00353.
https://rss.onlinelibrary.wiley.com/doi/...
) and the criterion of density predictive ordered - CPO (GELFAND & DEY, 1994GELFAND, A. & DEY , K. D. Bayesian Model Choice: Asymptotics and Exacts Calculations. Journal of the Royal Statistical Society: Series B, v.56, n.3, p.501-514, 1994. Available from: <Available from: https://rss.onlinelibrary.wiley.com/doi/abs/10.1111/j.2517-6161.1994.tb01996.x >. Accessed: Jun. 01, 2021. doi: 10.1111/j.2517-6161.1994.tb01996.x
https://rss.onlinelibrary.wiley.com/doi/...
) were used. As mentioned in ANDRADE FILHO et al. (2010ANDRADE FILHO, M. G. de et al. Bayesian approach to growth curve with objective priori. Revista Brasileira de Biometria, São Paulo, v.28, n.2, p.161-181, 2010. Available from: <Available from: http://jaguar.fcav.unesp.br/RME/fasciculos/v28/v28_n2/A10_Marinho.pdf >. Accessed: Feb. 01, 2021.
http://jaguar.fcav.unesp.br/RME/fascicul...
), lowest value for DIC and the highest for CPO indicated better fitted.

All the procedures necessary for this research, such as the generation of samples, summary of distributions, obtaining graphics, CPO and DIC, among others, were performed using the R software (R DEVELOPMENT CORE TEAM, 2020R DEVELOPMENT CORE TEAM. R: a language and environment for statistical computing. Vienna: R Foundation for Statistical Computing, 2020. Available from: <Available from: http://www.r-project.org >. Accessed: Jan. 20 2020.
http://www.r-project.org...
), mainly with packages nlme (PINHEIRO et al., 2019PINHEIRO J. et al. R Core Team. nlme: Linear and Nonlinear Mixed Effects Models. R package version 139, 2019. Available from: <Available from: https://CRAN.R-project.org/package=nlme >. Accessed: Feb. 01, 2021.
https://CRAN.R-project.org/package=nlme...
) and coda (PLUMMER et al., 2006PLUMMER M. et al. “CODA: Convergence Diagnosis and Output Analysis for MCMC.” R News, v.6, n.1, p.7-11, 2006. Available from: <Available from: https://journal.r-project.org/archive/ >. Accessed: Feb. 01, 2021.
https://journal.r-project.org/archive/...
).

RESULTS AND DISCUSSION:

Table 2 presents the estimates obtained for the mean and mode a posterior, as well as the HPD interval, for all parameters of the adjusted models, considering the non-irrigated regime (NI).

Table 2
Posterior estimates for the mean, mode and HPD interval of the parameters of the Logistic (L), Brody (B), von Bertalanffy (VB) and Richards (R) models, considering the data referring to the NI regime.

Analyzing table 2, it can be seen that, in general, the estimates for the mean and the mode a posterior of the parameters β and k for the Logistic, Brody and von Bertalanffy models are close, indicating that the a posterior marginal distribution of these parameters can be considered symmetrical. When analyzing the growth of the tomato plant fruit, MENDES et al. (2021MENDES, P. N. et al. Logistic Bayesian model in the study of the growth of tomatoes . Research, Society and Development, v.10, n.3, p.1-12, 2021. Available from: <Available from: https://rsdjournal.org/index.php/rsd/article/view/13198 >. Accessed: Apr. 05, 2021. doi: 10.33448/rsd-v10i3.13198
https://rsdjournal.org/index.php/rsd/art...
) also reported that the marginals of the parameters β and k, for the Logistic model are symmetrical.

However, when comparing the upper quantile of the HPD interval of the marginal distribution of the k parameter referring to the Brody model (qhpd), whose area below it is 97.5%, with the quantile of the Normal distribution (qn), considering k̅=0.00061 and sk = 0.00006 (not shown in the table) have: qhpd = 0.00071 and qn = 0.000727. Thus, as qn > qhpd it is concluded that the Normal has tails heavier than the marginal of parameter k.

Similarly, for parameter β of the Brody model, we have that qhpd = 0.97142 and qn = 0.974174 (considering β=0.9563 and sβ=0.00912 ). Therefore, as qn > qhpd it follows that Normal has tails heavier than the marginal of parameter β.

Regarding the estimates obtained for the parameters of the Richards model, there is a marked difference between the values of the mean and mode a posterior (Table 2). This indicate that the marginals of the parameters are asymmetric, and according to SAVIAN et al. (2009SAVIAN, T. V. et al. Bayesian analisys for ruminal degradability models. Ciência Rural, Santa Maria, v.39, n.7, p.2169-2177, out. 2009. Available from: <Available from: https://www.scielo.br/scielo.php?pid=S0103-84782009000700033&script=sci_abstract&tlng=pt >. Accessed: Feb. 01, 2021. doi: 10.1590/S0103-84782009000700033.
https://www.scielo.br/scielo.php?pid=S01...
) and SILVA et al. (2020SILVA, E. M. et al. Bayesian approach to the zinc extraction curve of soil with sewage sludge. Acta Scientiarum. Technology, v.42, n.1, p.1-9, 2020. Available from: <Available from: http://periodicos.uem.br/ojs/index.php/ActaSciTechnol/article/view/46893 >. Accessed: Feb. 01, 2021. 10.4025/actascitechnol.v42i1.46893
http://periodicos.uem.br/ojs/index.php/A...
), in future studies, this information can be considered in specification of prior distributions.

Still in relation to table 2, it can be seen that, in general, the HPD intervals obtained for the parameters were significant, affirming the efficiency of the Bayesian methodology in the estimation of parameters, reported in several studies, such as ANDRADE FILHO et al. (2010ANDRADE FILHO, M. G. de et al. Bayesian approach to growth curve with objective priori. Revista Brasileira de Biometria, São Paulo, v.28, n.2, p.161-181, 2010. Available from: <Available from: http://jaguar.fcav.unesp.br/RME/fasciculos/v28/v28_n2/A10_Marinho.pdf >. Accessed: Feb. 01, 2021.
http://jaguar.fcav.unesp.br/RME/fascicul...
), MARTINS FILHO et al. (2008MARTINS FILHO, S. et al. Bayesian approach in the growth curves of two cultivars of common bean. Ciência Rural, Santa Maria, v.38, n.6, p.1516-1521, 2008. Available from: <Available from: https://www.scielo.br/scielo.php?pid=S0103-84782008000600004&script=sci_abstract&tlng=pt >. Accessed: Feb. 01, 2021. doi: 10.1590/S0103-84782008000600004.
https://www.scielo.br/scielo.php?pid=S01...
) and SAVIAN et al. (2009SAVIAN, T. V. et al. Bayesian analisys for ruminal degradability models. Ciência Rural, Santa Maria, v.39, n.7, p.2169-2177, out. 2009. Available from: <Available from: https://www.scielo.br/scielo.php?pid=S0103-84782009000700033&script=sci_abstract&tlng=pt >. Accessed: Feb. 01, 2021. doi: 10.1590/S0103-84782009000700033.
https://www.scielo.br/scielo.php?pid=S01...
), among others.

Table 3 lists the estimates obtained for the mean and mode a posterior as well as the HPD interval, for all parameters of the adjusted models, considering the irrigated regime.

Table 3
Posterior estimates for the mean, mode and HPD interval of the parameters of the Logistic (L), Brody (B), von Bertalanffy (VB) and Richards (R) models, considering the data referring to the irrigated regime.

Based on table 3, it appears that, also in the irrigated regime, the HPD intervals were significant for all parameters. The interval obtained for parameter α, considering the adjustment of the Brody model, for example, indicates that the most likely maximum heights, for coffee trees at the adult phase, are between 210.6 cm and 234.8 cm.

Still in relation to Table 3, there was a marked difference between the estimates for the mean and mode a posterior of the parameter φ1 for the von Bertalanffy and Logistic models, suggesting that the marginals are asymmetric to the left. Furthermore, the HPD intervals for the φ parameter do not include zero (Table 2 and Table 3), emphasizing that the autocorrelation is significant in all studied scenarios. As highlighted by Silva et al. (2021SILVA, W. S. et al. Eucalyptus grandis x eucalyptus urophylla growth curve in different site classifications, considering residual autocorrelation. Revista Brasileira de Biometria, v.39, n.1, p.122-138, 2021. Available from: <Available from: https://biometria.ufla.br/index.php/BBJ/article/view/511 >. Accessed: Jun. 01, 2021. doi: 10.28951/rbb.v39i1.511
https://biometria.ufla.br/index.php/BBJ/...
), incorporating residual dependence may not necessarily improve the quality of the fit, but it is more consistent, since the assumption of independence is not satisfied in longitudinal data modeling.

The differences between the values obtained for the mean and mode a posterior for all parameters of the Richards model were also accentuated, indicating that the marginal distributions for its parameters are asymmetric. Such behavior was also verified for the marginal distributions of parameters of this model considering the NI regime (Table 2).

To compare the quality of fit provided by the four models analyzed, the Deviance information criterion (DIC) and the criterion of predictive ordered (CPO) density were calculated, the results are listed in table 4. As mentioned lowest value for DIC and the highest for CPO indicated better fitted.

Table 4
Estimates for the DIC and CPO selection criteria, considering data referring to the non-irrigated (NI) and irrigated regimes and the Logistic (L), Brody (B), von Bertalanffy (VB) and Richards (R) models.

Considering the NI regime, the estimates obtained for both adopted criteria indicate that, among the models evaluated, the Logistic model best describes the growth in height of coffee trees, as this model presented the lowest value for DIC and the highest for CPO (Table 4). The marginal distributions of each parameter of the Logistic model for the NI regime are illustrated in figure 1.

Figure 1
Sample histograms of the marginal distributions of parameters of the Logistic model considering first order residual autocorrelation and the NI regime.

Regarding the irrigated regime, the CPO and DIC estimates indicated that the Brody model, among the evaluated models, is the one that best represents the growth curve in height of the irrigated coffee tree (Table 4). The marginal distributions of each parameter of the Brody model for the irrigated regime are illustrated in figure 2.

Figure 2
Sample histograms of the marginal distributions of parameters of the Brody model considering data referring to the irrigated regime.

The growth of the NI and irrigated coffee trees followed different growth patterns. As the Logistic model was more suitable to describe the NI coffee tree height, sigmoidal growth was observed with maximum growth rate at 726 days after planting (estimate of parameter β of the Logistic model in Table 2). For the irrigated coffee tree, the Brody model was more suitable, indicating that the growth started at the maximum growth rate decreasing over time (SANTOS et al., 2019SANTOS, H. B. et al. Application of non-linear mixed models for modelling the quail growth curve for meat and laying. The Journal of Agricultural Science, v.156, n.10, p.1216-1221, 2019. Available from: <Available from: https://www.cambridge.org/core/journals/journal-of-agricultural-science/article/abs/application-of-nonlinear-mixed-models-for-modelling-the-quail-growth-curve-for-meat-and-laying/88C026C804B4297449AB4539D2AADAA3 >. Accessed: Feb. 01, 2021. doi: 10.1017/S0021859619000169
https://www.cambridge.org/core/journals/...
). Moreover, the HPD for the α parameter of the Logistic model for NI coffee was 171cm to 186 cm (Table 2) and for the Brody model for irrigated coffee tree it was 210 cm to 234 cm (Table 3), thus indicating that the irrigated coffee tree had a maximum height greater than the non-irrigated plant, as there was no intersection in the intervals. Importantly, one of the great advantages of working with the Bayesian methodology is to obtain the credibility interval (HPD), which is the interval with (1 - α)% probability that contains the most plausible values for the parameter, which allows direct comparison between intervals (GUEDES et al., 2014GUEDES, T. A. et al. Nonlinear models applied to seed germination of Rhipsalis cereuscula Haw (Cactaceae). Acta Scientiarum. Technology, v.36, n.4, p.651-656, 2014. Available from: <Available from: http://periodicos.uem.br/ojs/index.php/ActaSciTechnol/article/view/21192 >. Accessed: Feb. 01, 2021. doi: 10.4025/actascitechnol.v36i4.21192.
http://periodicos.uem.br/ojs/index.php/A...
; SILVA et al., 2020SILVA, E. M. et al. Bayesian approach to the zinc extraction curve of soil with sewage sludge. Acta Scientiarum. Technology, v.42, n.1, p.1-9, 2020. Available from: <Available from: http://periodicos.uem.br/ojs/index.php/ActaSciTechnol/article/view/46893 >. Accessed: Feb. 01, 2021. 10.4025/actascitechnol.v42i1.46893
http://periodicos.uem.br/ojs/index.php/A...
).

The growth pattern, as well as the values observed for the mean height of non-irrigated coffee trees and the graphical representation of the Bayesian adjustment of the Logistic model with AR(1), as well as the values observed for the irrigated coffee trees and the Bayesian adjustment of the Brody model are shown in figure 3.

Figure 3
Mean height of coffee trees and Bayesian fitting of the Logistic with AR(1) and Brody models considering the NI and irrigated regimes, respectively.

CONCLUSION:

For non-irrigated plants, the Bayesian information criterion and the criterion of predictive ordered density indicated that, among the evaluated models, the Logistic model is the one that best described the coffee plant height growth over time. For irrigated plants, these same criteria indicated the Brody model. Thus, the growth of the non-irrigated and irrigated coffee trees followed different growth patterns, the NI coffee height showed sigmoidal growth with maximum growth rate at 726 days after planting and the irrigated coffee tree starts its development with high growth rates that gradually decrease over time.

ACKNOWLEDGEMENTS

“And was financed in part by the Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq).”

REFERENCES

  • CR-2021-0275.R1

Edited by

Editors:

Leandro Souza da Silva(0000-0002-1636-6643) Alberto Cargnelutti Filho(0000-0002-8608-9960)

Publication Dates

  • Publication in this collection
    14 Mar 2022
  • Date of issue
    2022

History

  • Received
    08 Apr 2021
  • Accepted
    10 Oct 2021
  • Reviewed
    10 Dec 2021
Universidade Federal de Santa Maria Universidade Federal de Santa Maria, Centro de Ciências Rurais , 97105-900 Santa Maria RS Brazil , Tel.: +55 55 3220-8698 , Fax: +55 55 3220-8695 - Santa Maria - RS - Brazil
E-mail: cienciarural@mail.ufsm.br