# Source paper archive - Source title: A hierarchical Bayesian model for dynamic phenotypes in multi-environment trials with genotype by environment interaction - Author(s): Shota Gugushvili, George A. K. van Voorn, Martin P. Boer, Daniela Bustos-Korts, Sandra Huynh Truong, Nicholas A. Friedenberg, Ryan McCormick, Scott Chapman, Carlos D. Messina, Fred A. van Eeuwijk - Source landing URL: https://www.biorxiv.org/content/10.64898/2026.11.05.764527v1 - PDF URL: https://www.biorxiv.org/content/10.64898/2026.11.05.764527v1.full.pdf - PDF SHA-256: d6cd6cdac6ae4dc41edbf6028fd9bd3c0ce6fc7debf8126e1d277446606cacee - Converter: pymupdf4llm 0.0.27 - Limits: Automated PDF text extraction; reading order, equations, tables and formatting may be incomplete or incorrect. Images are omitted, and no OCR was performed. Markup and links are neutralized for safe downloading. Consult the original PDF for the authoritative version. --- [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.05.764527; this version posted November 6, 2026. The copyright holder for this](https://doi.org/10.64898/2026.11.05.764527) preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. # 1 A hierarchical Bayesian model for dynamic phenotypes in 2 multi-environment trials with genotype by environment interaction 3 Shota Gugushvili [1,*], George A. K. van Voorn [1], Martin P. Boer [1], Daniela Bustos-Korts [1,2], 4 Sandra Huynh Truong [3], Nicholas A. Friedenberg [3], Ryan McCormick [3], Scott Chapman [4], 5 Carlos D. Messina [3,5], and Fred A. van Eeuwijk [1,*] 6 1 Mathematical and Statistical Methods group – Biometris, Wageningen University and Research, Wageningen, 7 The Netherlands 8 2 Institute of Plant Production and Protection, Faculty of Agricultural Sciences, Universidad Austral de Chile, 9 Valdivia, Chile 3 10 Research & Development, Corteva Agriscience, Johnston, IA, United States 4 11 School of Agriculture and Food Sustainability, Faculty of Science, The University of Queensland, Gatton, QLD, 12 Australia 13 5 Department of Horticultural Sciences, University of Florida, Gainesville, FL, United States 14 Corresponding authors: [shota.gugushvili@wur.nl,](mailto:shota.gugushvili@wur.nl) [fred.vaneeuwijk@wur.nl](mailto:fred.vaneeuwijk@wur.nl) 15 **Abstract** 16 New phenotyping tools produce longitudinal phenotypes that describe plant growth and devel17 opment. Estimates for the parameters that describe the dynamics of such longitudinal phenotypes 18 can be incorporated as predictors in models for static target traits like yield and biomass. Within a 19 hierarchical Bayesian context, we present a differential equation (DE) model for log-logistic growth in 20 longitudinal traits for sets of genotypes across multiple environments. The DE parameters coincide 21 with a midpoint, a slope, and an asymptote at the end of growth and development, i.e., final biomass. 22 Each DE parameter in our hierarchical model results from contributions of genotype, environment, 23 and _G × E_ . This formulation allows one to assess and quantify the genotype-specificity of the DE 24 parameters. When none of the DE parameters in our model requires a _G × E_ term, the dynamic 25 parameters characterizing the particular trait will be useful features for inclusion in prediction models 26 for yield or biomass. We fitted a DE model to APSIM-simulated longitudinal wheat biomass data 27 and quantified the contribution of genotype, environment and _G × E_ to the phenotypic variation in 28 the DE parameters. Only the upper asymptote showed clear _G × E_ . Estimates for genotypic and 29 _G × E_ parameters within the DE parameters were used to characterize genotypes. Estimates for 30 environmental and _G × E_ effects formed inputs for environmental characterization. 31 **Keywords:** crop growth model, genotype by environment interaction, phenomics, phenotyping, 32 process-based model, yield prediction 1 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.05.764527; this version posted November 6, 2026. The copyright holder for this](https://doi.org/10.64898/2026.11.05.764527) preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. ## 33 1 Introduction 34 New phenotyping devices produce longitudinal data on many aspects of plant growth and development for 35 large collections of genotypes (Gu´edon et al., 2022). Depending on the device, these data have moderate to 36 high time resolution (Araus et al., 2018, 2022; Tardieu et al., 2017). The study of longitudinal phenotypic 37 traits is of interest in itself when these traits can be considered to represent yield components, i.e. they 38 bear a physiological relation to yield and reflect adaptation of the genotype to the environment. Knowledge 39 of the dynamics of yield components adds to biological insights in yield-producing processes and will 40 improve yield prediction models. Such improvements can cover both within-trial yield forecasts as well as 41 between-trial yield predictions where the longitudinal phenotyping information can be combined with 42 genomic information. 43 Beyond a small subset of longitudinal traits that can be considered to be directly or closely related to 44 yield, most longitudinal traits that are measured automatically or semi-automatically by proximate and 45 remote sensors initially lack a direct connection with plant physiological processes and structures. Such 46 traits are called secondary traits, and although they lack clear biological meaning, it can be worthwhile 47 to model their dynamics and estimate the corresponding parameters because aspects of their dynamical 48 behavior may improve prediction models for yield. When secondary traits or functions of secondary traits 49 contribute to the prediction of yield, then it is plausible that they represent proxies of physiologically 50 important traits, that is, yield components. Various articles offer guidelines on how to decompose yield 51 into components and which traits to phenotype as proxies to those components (Araus et al., 2018; 52 Tardieu et al., 2017; Bouidghaghen et al., 2023). These yield components characteristically reflect resource 53 capture and conversion into biomass and aspects of source-sink relations that lead to the prediction of the 54 fraction of biomass translated into yield (harvest index). A decomposition of yield into component traits 55 or processes immediately suggests yield prediction models that integrate information from multiple yield 56 components as measured by proxies. The yield components act as genotype-specific inputs. The prediction 57 model consists of a function that is able to convert the yield components together with some environmental 58 inputs into yield as output. Prediction models differ with respect to the yield components to be included 59 as static or dynamic quantities. The static yield components represent physiological constants that may be 60 observed in single well-conditioned experiments, while the dynamic yield components require longitudinal 61 observations. The environmental inputs are dynamical, but may differ at the aggregation level that is 62 incorporated. For example, temperatures may be chosen to be included as a daily or weekly mean. 63 Bouidghaghen et al. (2023) make a case for the breakdown of yield into a number of subprocesses and 64 structures, each of which can be measured by high-throughput robotized phenotyping devices at platforms. 65 They emphasize that before inserting traits measured on platforms into prediction models for yield in the 66 field, these traits require scaling, transformation, and conditioning on environmental variables such as, for 67 example, thermal time, soil water potential, and vapor pressure deficit. Although Bouidghaghen et al. 68 (2023) illustrate how to decompose yield and predict yield components from phenotyping and genomics, 69 they do not go all the way in formulating and comparing prediction models for yield from genotype-specific 70 and environmental inputs. The same approach to yield prediction via yield decomposition and component 71 approximation drives the work of other researchers that aims at the integration of longitudinal phenotyping 72 information and ultimately genomics into yield prediction models. Examples of this approach can be 73 found in Roth et al. (2021, 2022, 2023) and P´erez-Valencia et al. (2022, 2024). All of these papers focus 74 on one or a few longitudinal traits whose dynamics are modeled to extract summary statistics for eventual 75 prediction of yield in future environments. Again, this latter aspect is not explicitly elaborated. 76 The papers mentioned in the above paragraphs all depart from dynamical models for yield (or yield 77 components) that are based on physiological insights and help to define underlying yield components 78 that can either be measured by sensors in the field or on a platform or predicted from marker profiles 79 in an attempt to integrate phenotyping with genomics in the prediction of yield. The dynamical model 80 can be a fully developed ecophysiological model, such as in Technow et al. (2015) and Messina et al. 81 (2018), or more data driven as in the case of Roth et al. (2021, 2022) and P´erez-Valencia et al. (2022, 82 2024). In all cases, an important consideration in the use of longitudinal phenotypic information for 83 yield prediction is that the dynamics of the longitudinal trait itself, say the defining parameters, are not 84 subject to strong genotype by environment interaction, _G × E_, themselves (Poudel et al., 2023; Roth 85 et al., 2023). The idea here is that parameters describing longitudinal yield components will enter as 86 genotype-specific static inputs into yield generating functions together with dynamical environmental 87 inputs. When genotype-specific inputs do not show _G × E_, or when they show a simple type of _G × E_ 88 that can be scaled, transformed, or conditioned on an environmental covariable (Bouidghaghen et al., 89 2023), they can be measured or estimated in a single platform or field experiment. _G × E_ for yield then 2 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.05.764527; this version posted November 6, 2026. The copyright holder for this](https://doi.org/10.64898/2026.11.05.764527) preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 90 needs to emerge from the yield generating function. 91 The objective of this paper is to investigate a class of statistical models for analyzing longitudinal 92 trait data for a population of genotypes observed by one or more phenotyping devices in a series of 93 environments. The longitudinal traits may have been chosen after a prior yield decomposition and the 94 expectation is that the parameters describing the longitudinal behavior of the component trait will exhibit 95 relatively little _G × E_ and that they will ultimately contribute to a better yield prediction. For improved 96 yield prediction, we may later need to combine dynamic information on multiple yield components and 97 environmental inputs, but on the road to a robust predictive model, we first need to develop adequate 98 models for describing individual longitudinal traits. We aim for a modeling framework for our individual 99 longitudinal traits that is flexible in the kind of functions and parameterizations that can be chosen to 100 describe the dynamics of the trait. Simultaneously, across the entire set of genotypes and environments, 101 we would like to test and quantify the amount of _G × E_ in each of our dynamic parameters, using the 102 same inferential scheme that was used to estimate these parameters. 103 We will elaborate here the modeling concepts for plant growth that were introduced in Van Voorn et al. 104 (2023). In that paper, simple differential equation (DE) models with few state variables and parameters 105 were presented that were derived from first principles. Subsequently, these models were fitted per genotype 106 and per trial, that is, the estimation process was done separately for each combination of genotype and 107 environment without borrowing information between genotypes or environments, although it can be 108 assumed that a population of genotypes will show similar growth patterns over time just as environments 109 with comparable stress will do. Fitting dynamic models for single genotype by environment combinations 110 is feasible when the time resolution for observations is high, say, daily to a few times per week (Roth 111 et al., 2021, 2022; Shao et al., 2026). However, such a resolution is not realistic for most field situations. 112 Longitudinal trait data produced by drones in outside-field experiments typically consist of series of 10 to 113 15 time points in a growing season. For such medium-resolution temporal data, borrowing information 114 between genotypes in the fitting process will be beneficial, and very likely the same is true for borrowing 115 information between environments. With medium-resolution secondary phenotypes in mind, we present 116 here a hierarchical Bayesian framework for the simultaneous fitting of dynamic models on longitudinal 117 trait data for populations of genotypes across multiple environments. This framework allows borrowing 118 of information between genotypes and environments and has as further advantages the use of shrinkage 119 estimation of the dynamical parameters as well as appropriate propagation of uncertainty on parameter 120 estimates and phenotypic predicted values. As the most distinguishing feature, our framework expresses 121 dynamic parameters as functions of genotypic, environmental, and _G × E_ terms in a fashion reminiscent of 122 generalized linear mixed models (GLMMs). The decomposition of the phenotypic variation of the dynamic 123 parameters into components attributable to genotype, environment, and _G × E_ provides information on 124 the genotype specificity of these parameters and the genotypic and environmental drivers of dynamic 125 behavior. An example analysis on APSIM simulated wheat ( _Triticum_ _aestivum_ L.) biomass data, the 126 generation of which is described in Bustos-Korts et al. (2019a,b), illustrates the machinery of our Bayesian 127 hierarchical approach together with the interpretation of the estimates of the genotypic and environmental 128 contributions to the DE parameters. ## 129 2 Materials and methods 130 In this section, we present key aspects of our approach. We first give a rough and intuitive sketch of our 131 approach and philosophy. We then recall a phenomenological model from Van Voorn et al. (2023) and 132 show how a modification of it can be embedded within a statistical model. The resulting model will be 133 fitted via a Bayesian approach. An APSIM simulated wheat dataset is used to illustrate and investigate 134 the practical performance of our method. 135 **2.1** **A** **sketch** **of** **our** **hierarchical** **approach** 136 In a generic notation, we can write the signal part of the dynamics of the _s_ _[th]_ secondary trait for a genotype, 137 _i_, in an environment, _j_, as: d _x_ _ij,s_ ( _t_ ) = _f_ ( _x_ _ij,s_ ( _t_ ) _, θ_ _ij_ _, z_ _j_ ( _t_ )) _._ d _t_ 138 Here _x_ _ij,s_ ( _t_ ) represents, for genotype _i_ in environment _j_, an error free time-dynamic secondary trait, _s_, 139 such as biomass. The variable _z_ _j_ ( _t_ ) describes a vector of time-dynamical exogenous inputs, such as one or 140 more soil properties, minimum, mean and maximum temperatures, and precipitation. The DE or process 3 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.05.764527; this version posted November 6, 2026. The copyright holder for this](https://doi.org/10.64898/2026.11.05.764527) preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. Table 1: Explanations of the symbols and abbreviations used throughout this paper. **Symbol** **Meaning** CGM (abbr.) Crop Growth Model DE (abbr.) Differential Equation GLMM (abbr.) Generalized Linear Mixed Model _G × E_ (abbr.) Genotype-by-Environment Interaction WGP (abbr.) Whole Genome Prediction _t_ Continuous time variable _i_ Index for genotype _j_ Index for trial (location _×_ year) _k_ Index for time (discrete) _p_ Superscript for DE parameter _s_ Subscript to indicate secondary trait _x_ _ij,s_ ( _t_ ) Secondary trait state variable _y_ _ijk_ Observation at time _t_ _k_ _z_ _j_ ( _t_ ) Exogenous (environmental) input _σ_ Standard deviation _θ_ _ij_ Parameter (generic) _µ_ Intercept for DE parameter _α_ _i_ Genotypic effect on DE parameter _β_ _j_ Trial effect on DE parameter ( _α · β_ ) _ij_ Interaction effect on DE parameter _M_ ( _t_ ) Total biomass (time-dependent) _r_ Crop growth rate _M_ _[max]_ Maximum attainable biomass _t_ _[mid]_ Midpoint _∼_ Symbol for ‘distributed according to’ _N_ ( _µ, σ_ ) Normal distribution with mean _µ_ and standard deviation _σ_ Exponential( _λ_ ) Exponential distribution with rate parameter _λ_ HalfNormal( _µ, σ_ ) Half-Normal distribution with parameters _µ_ _σ_ 4 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.05.764527; this version posted November 6, 2026. The copyright holder for this](https://doi.org/10.64898/2026.11.05.764527) preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 141 parameters are contained in the vector _θ_ _ij_ . From this ingredient we build a nonlinear mixed model, in 142 which each DE parameter is related to a linear mixed predictor through a log link, in the spirit of a 143 GLMM. In particular, the hierarchical system we propose can be thought of as containing an extra level 144 of modeling in defining GLMMs for each of the DE parameters: _g_ ( _θ_ _ij_ ) = _Xγ_, where the link function 145 _g_ ( _·_ ) works on the expectation of _θ_ _ij_ and is linear in the regression parameters _γ_, with design matrix _X_ . 146 The GLMMs for the DE parameters will include a general intercept, a genotypic main effect term, an 147 environmental main effect term, and a _G × E_ term. This formulation allows testing explicitly for the 148 inclusion of _G × E_ terms in the GLMMs for the individual DE parameters. Subsets of genotypes and 149 environments for whose dynamics no _G_ _×_ _E_ terms are needed act in what we can call a genetically additive 150 way, i.e., it suffices to estimate the differences between genotypes in a single environment because the same 151 genotypic differences will occur in any environment in which the genotypes will be grown, i.e. predictions 152 are valid for any environment. To obtain a genomic prediction model, the genotypic main effects can 153 be replaced or approximated by penalized regressions on marker profiles. It is important to note that 154 additivity of genotype and environment for the DE parameters does not exclude _G × E_ at the level of the 155 dynamic trait itself because the non-linearity of the DE system can still create _G × E_ as an emerging 156 property on top of DE parameters that behave additively with respect to genotype and environment. 157 **2.2** **Phenomenological** **model** 158 Van Voorn et al. (2023) discuss four models that reflect some essential characteristics of crop growth. 159 These models are of low complexity, while being time-dynamical, as well as phenomenological in nature, 160 and offer a biological interpretation of their parameters. According to Van Voorn et al. (2023), genotypic 161 variation can be conceptualized as genotype-related differences in parameter values, while dependence 162 on environmental or management inputs can be explicitly introduced by making the model parameters 163 functions of such inputs. 164 In this paper, we focus on a logistic model for biomass in wheat. This is the most basic model in 165 Van Voorn et al. (2023), which we will combine with a hierarchical Bayesian modeling approach. The 166 resulting statistical model is flexible enough to reach non-trivial inferential conclusions. In addition, 167 the logistic model serves as a useful baseline against which other approaches can be benchmarked and 168 compared. Furthermore, there is room to explore common numerical challenges in nonlinear regression 169 and fitting of ordinary differential equation (ODE) models (Seber and Wild, 1989; Ramsay and Hooker, 170 2017), as well as to examine the scalability of the method to large data. 171 **2.2.1** **Logistic** **model** 172 In the notation of Van Voorn et al. (2023), the logistic model is given by the ODE 1 _−_ _[M]_ [(] _[t]_ [)] _M_ _[max]_ )︃ d _M_ ( _t_ ) = _rM_ ( _t_ ) d _t_ (︃ 173 with an explicit solution _M_ _[max]_ _M_ ( _t_ ) = 1 + ( _,_ _M_ (0) = _M_ 0 _,_ (1) _._ (2) exp( _−rt_ ) _M_ _[max]_ ) _−_ 1 _M_ 0 174 Here _M_ ( _t_ ) is the total biomass at time _t_, the positive parameter _r_ is an intrinsic growth rate, while the 175 positive parameter _M_ _[max]_ is an intrinsic growth limitation or asymptote. Depending on the context, the 176 initial condition _M_ 0 can be interpreted either as the seed weight, or as the biomass at the starting day of 177 measurement. 178 In the literature there exist several alternative parameterizations of the logistic model (1) – (2), see, e.g., 179 Seber and Wild (1989) and Pinheiro and Bates (2000). In one parameterization, the three parameters of 180 the model are the intrinsic growth rate _r_, the intrinsic growth limitation or asymptote _M_ _[max]_, and the 181 midpoint. _−_ 1 )︃ _t_ _[mid]_ = [1] _r_ [log] _M_ _max_ (︃ _M_ 0 _,_ 182 where the biomass _M_ ( _t_ _[mid]_ ) = _M_ _[max]_ _/_ 2 _._ Here we implicitly assume that _M_ 0 _< M_ _[max]_ . The function itself 183 is then given by _M_ _[max]_ _M_ ( _t_ ) = 1 + exp( _−r_ ( _t −_ _t_ _[mid]_ )) _[.]_ (3) 184 This parameterization will be the basis for our subsequent arguments. 5 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.05.764527; this version posted November 6, 2026. The copyright holder for this](https://doi.org/10.64898/2026.11.05.764527) preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 185 **2.2.2** **Log-logistic** **model** 186 Successful fitting of nonlinear regression models to the data is strongly dependent on the chosen parameter187 ization of the function to be fitted (Seber and Wild, 1989, pp. 673–674). For example, a reparameterization 188 of a nonlinear model could remedy parameter identifiability (Seber and Wild, 1989, p. 103). A suitable 189 reparameterization of a model could likewise improve the convergence of an iterative non-linear optimiza190 tion algorithm (Seber and Wild, 1989, pp. 95–97), or, in the context of the present paper, convergence of 191 an MCMC sampler. Ross (1990) contains further discussion on these topics. 192 A modified version of the logistic model (1) – (2), termed the log-logistic model, turned out to be 193 numerically easier and faster to fit in the hierarchical modeling context than the basic logistic model. This 194 appears to be caused by algorithmic peculiarities of our MCMC sampler. In general, both the logistic and 195 log-logistic curves can be expected to give visually quite similar fits. 196 The log-logistic growth function closely resembles the logistic model (3) and is given by ︁ _M_ _[max]_ _M_ ( _t_ ) = ︁ )︁ _−r_ _[.]_ (4) 1 + ( _t_ _[mid]_ _t_ ︁ ︁ 197 The three parameters of the log-logistic model have an intuitive interpretation similar to the parameters of 198 the logistic model: _M_ _[max]_ is the intrinsic growth limitation. Next, _t_ _[mid]_ is the time point where response 199 _M_ ( _t_ _[mid]_ ) = _M_ _[max]_ _/_ 2. Finally, the parameter _r_ is the scaled slope at the midpoint. See Ritz et al. (2020) 200 for additional details. 201 The log-logistic model is by far the most common dose response model (Ritz et al., 2020, p. 178). 202 Its name reflects the fact that the function is obtained by rescaling the log-logistic distribution, which 203 in economics is known as the Fisk distribution. The log-logistic distribution also arises in the survival 204 analysis context (Bennett, 1983). Its mathematical properties make it suitable for crop growth modeling 205 as well, and in fact the model has been used in the agricultural setting: some references are listed in Focke 206 et al. (2017). 207 The differential equation for the log-logistic growth model away from the time origin _t_ = 0 can be 208 shown to be d _M_ ( _t_ ) _[r]_ _[M]_ [(] _[t]_ [)] ︁ )︃ ︁ = _[r]_ d _t_ _t_ ︁ _t_ _[M]_ [(] _[t]_ [)] ︁ (︃ ︁ 1 _−_ _[M]_ [(] _[t]_ [)] _M_ _[max]_ ︁ _._ (5) ︁ 209 The initial condition is _M_ ( _t_ _[mid]_ ) = _M_ _[max]_ _/_ 2 _,_ and the equation is thus understood to operate forward and 210 backward in time. When fitting the log-logistic model to data, we will rely on the explicit form (4) rather 211 than the differential equation formulation (5). See Section 2.5 ahead. 212 **2.3** **Data** **description** 213 To illustrate our method, we chose the same simulated longitudinal wheat data as in Van Voorn et al. (2023). 214 The crop growth simulation platform was APSIM-Wheat (Keating et al., 2003; Holzworth et al., 2018).The 215 data consist of biweekly (i.e. successive measurements spaced two weeks apart) biomass observations over 216 the growing season for an Australian diversity panel with originally 199 wheat genotypes (Zheng et al., 217 2013) observed at four locations over 31 growing seasons (years); see Bustos-Korts et al. (2019a) and 218 Bustos-Korts et al. (2019b). The biweekly observations were chosen to reflect typical medium-resolution 219 temporal data that are similar to drone imagery observations for biomass. The complete data set is 220 available online (Bustos-Korts, 2023). For computational convenience, we drew a genetically representative 221 subset of 48 genotypes from the entire set of 199 genotypes (see the code referred to in the Software 222 section). This subsetting of genotypes allowed us to do all calculations described in this paper on a laptop 223 (MacBook Air, 2020, with Apple M1 chip with 8 cores, 4 for performance and 4 for efficiency, and 8 GB 224 RAM). The four locations chosen (Emerald, Narrabri, Yanco, Merredin) corresponded to environmentally 225 contrasting wheat growing areas in the Australian wheat belt. The environments in this paper are a 226 subset of the environments described in Casadebaig et al. (2016) and Chenu et al. (2013). 227 The first 14 days of biomass data after sowing were discarded, as these biomasses were (near) zero, 228 which presented numerical problems for fitting. The biomass given in kg ha _[−]_ [1] was converted into tons 229 ha _[−]_ [1], and the time in days was re-scaled by a factor of 10, so that a new time unit corresponded to 10 230 days. These transformations bring both variables to the same scale and can be expected to help with 231 non-linear fitting, for instance because the population level parameters are then of the same order of 232 magnitude (Seber and Wild, 1989, p. 674). 233 At this stage of research, we preferred simulated data over field trial data. This choice allows for an 234 objective assessment of our approach against the ground truth. On the environmental side, the data were 6 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.05.764527; this version posted November 6, 2026. The copyright holder for this](https://doi.org/10.64898/2026.11.05.764527) preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 235 qualitatively characterized by the type of limitation (e.g. temperature or water) and quantitatively by 236 the level of environmental stress. On the genetic side, the data showed qualitative differences induced by 237 genotypes. Bustos-Korts et al. (2019a) and Bustos-Korts et al. (2019b) chose genotypic specifications of 238 crop growth parameters that mimic realistic genetic variation for Australian environmental conditions. 239 **2.4** **Bayesian** **approach** 240 Hierarchical models, of which our statistical model is a representative, lend themselves naturally to a 241 Bayesian approach. Moreover, in our experience with hierarchical models, modern Bayesian software such 242 as **Stan** is less prone to numerical failure than the off-the-shelf nonlinear least squares methods in **R** that 243 were used by Van Voorn et al. (2023). See Comets et al. (2017), Driver and Voelkle (2018) and Plan et al. 244 (2012) for a detailed discussion on algorithmic convergence issues in this context. 245 A data generating mechanism in a Bayesian model is described through the likelihood function. In 246 line with typical assumptions of additive and normally distributed measurement errors, we assumed that 247 the likelihood function was Gaussian and centered on values _M_ ( _t_ _k_ ) at observation points _t_ _k_, and had a 248 constant standard deviation _σ._ 249 Another mandatory ingredient of Bayesian approach is a prior distribution on the model parameters. 250 This is a probability distribution that reflects one’s prior knowledge on model parameters. It can also be 251 viewed as a penalty or regularization term in the estimation procedure. Transtrum et al. (2011) explain 252 the usefulness of priors as penalty terms in the context of non-linear least squares estimation. Our exact 253 prior specifications will be given later on. 254 Except for simple statistical models, the posterior distribution is rarely available in analytic form. 255 To address this issue, over the years powerful Markov Chain Monte Carlo (MCMC) methods have been 256 developed to approximate posterior distribution through random sampling; see Brooks et al. (2011) for 257 an overview. In general, MCMC methods are regarded as the gold standard for approximate posterior 258 inference in that, in a long run, they provide asymptotically exact posterior samples. Their accuracy is 259 offset by constraints of scalability and computational feasibility. Other options to approximate the posterior 260 include variational inference (Blei et al., 2017), sequential Monte Carlo (Chopin and Papaspiliopoulos, 261 2020), Approximate Bayesian Computation (ABC) (Sisson et al., 2019), and the integrated nested Laplace 262 approximation (INLA) (G´omez-Rubio, 2020). 263 **2.5** **Bayesian** **hierarchical** **model** 264 The log-logistic model (4) was used to develop a hierarchical Bayesian model (Gelman and Hill, 2006; 265 Gelman et al., 2014) for the APSIM wheat data from Section 2.3. The underlying idea of hierarchical 266 modeling is that even though genotypes differ from each other, they have common features, and one 267 can attempt to exploit this when fitting a model to the data. A similar idea applies to environmental 268 variations over multiple years. As a result, data are used more efficiently, which, in principle, should 269 lead to more accurate parameter estimates and a better quantification of the remaining uncertainty 270 in inferential conclusions. This becomes particularly important in the real data setting where biomass 271 measurements are obtained through relatively infrequent drone flights; see Chen et al. (2022), Smith et al. 272 (2021), Panday et al. (2020a), Panday et al. (2020b) for an overview of longitudinal drone imagery. 273 In a nutshell, we assume that dynamics of crop growth in each trial and for each genotype can be 274 described by a log-logistic function with an own set of parameters, whereas the remaining variation is 275 attributed to statistical noise. The three parameters of the log-logistic function are themselves driven 276 by genetic and environmental information. In the next step, we assume that parameters across trials 277 and genotypes are linked through an overarching distribution. Inference is to be conducted on this 278 latter distribution. The resulting statistical model is dynamic, but embedded within the mixed modeling 279 framework at the parameter level (Pinheiro and Bates, 2000; Lee, 2022). 280 For the model description below, we provide a Kruschke-style diagram, cf. Kruschke (2014), in Figure 281 1. We use the following notation throughout: _y_ _ijk_ is the biomass measurement at time _t_ _k_ for genotype _i_ 282 and trial _j._ Next, _M_ _ijk_ = _M_ ( _t_ _k_ ; _θ_ _ij_ ) is the value of the log-logistic function at time _t_ _k_ for genotype _i_ and 283 trial _j._ The parameter _θ_ _ij_ is a three-dimensional vector, its coordinates being logarithms of the maximum 284 or right asymptotic biomass _M_ _ij_ _[max]_, midpoint time _t_ _[mid]_ _ij_ and growth rate _r_ _ij_, respectively (for economy of 285 notation, here we have redefined parameters _θ_ _ij_ from Section 2.1 to their logarithms). Finally, _σ_ is the 286 measurement noise level. 287 As the construction is the same for each of the parameters of the log-logistic function, we focus on the 7 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.05.764527; this version posted November 6, 2026. The copyright holder for this](https://doi.org/10.64898/2026.11.05.764527) preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 288 right asymptote _M_ _[max]_ in our presentation. Our key assumption is the decomposition _ij_ _[max]_ = exp( _µ_ _[M]_ + _α_ _i_ _[M]_ _j_ _[M]_ + ( _α · β_ ) _[M]_ _ij_ _ij_ [)] _[.]_ (6) _M_ _[max]_ _i_ _[M]_ + _β_ _[M]_ 289 The exponentiation is to ensure positivity of the parameter, whereas working on the logarithmic scale 290 enables one to build into it a two-way analysis of variance or ANOVA model structure _θ_ _[M]_ _ij_ _[M]_ [=] _[ µ]_ _[M]_ [+] _[ α]_ _i_ _[M]_ _i_ _[M]_ + _β_ _[M]_ _j_ _[M]_ + ( _α · β_ ) _[M]_ _ij_ _ij_ 291 with, obviously, _θ_ _ij_ _[M]_ [=] [ log] [(] _[M]_ _[ max]_ _ij_ ). Next, _µ_ _[M]_ can be thought of as the intercept common to all genotypes 292 and trials. The term _α_ _i_ _[M]_ is contribution of the genotype _i_ (on the logarithmic scale) to the DE parameter 293 _M_ _ij_ _[max]_ . Similarly, _β_ _j_ _[M]_ models the effect of specific environmental and management conditions in the trial _j_ 294 on the DE parameter _M_ _ij_ _[max]_ . Finally, ( _α · β_ ) _[M]_ _ij_ [is] [the] [interaction] [term] [for] [genotype] _[i]_ [and] [trial] _[j.]_ [We] [treat] 295 _α_ _i_ _[M]_ [’] [s,] _[β]_ _j_ _[M]_ [’] [s] [and] [(] _[α][ ·][ β]_ [)] _[M]_ _ij_ [’] [s] [as] [random] [effects;] [see,] [for] [example,] [Pinheiro] [and] [Bates] [(][2000][).] [The] [linear] 296 model formulation for the logarithm of _M_ _ij_ _[max]_ is an example of a GLMM (Searle et al., 2006; Faraway, 297 2016) for a DE parameter embedded in our hierarchical Bayesian framework for multi-genotype multi298 environment dynamic data. This GLMM-like embedding offers a very powerful and versatile modeling 299 framework with options for quantification of genotypic and environmental contributions to dynamic trait 300 behavior. 301 The interpretability of the parameters of the log-logistic function can be directly leveraged to set priors 302 in the Bayesian model. The intercept _µ_ _[M]_ is assigned the prior with, obviously, _θ_ _[M]_ _ij_ _[M]_ [=] [ log] [(] _[M]_ _[ max]_ _ij_ _ij_ _[max]_ . Similarly, _β_ _j_ _[M]_ _ij_ _[max]_ . Finally, ( _α · β_ ) _[M]_ _ij_ _i_ _[M]_ [’] [s,] _[β]_ _[M]_ _j_ _[M]_ [’] [s] [and] [(] _[α][ ·][ β]_ [)] _[M]_ _ij_ _µ_ _[M]_ _∼_ _N_ (2 _,_ 2) _._ 303 This prior is weakly informative: the distribution covers a likely range of the logarithm of the biomass 304 measured at the end of the growing season for the APSIM data, while still allowing for a significant 305 variation. 306 The genotype random effects _α_ _i_ _[M]_ are assigned a prior through a hierarchical construction: _i_ _[M]_ _∼_ _N_ (0 _, σ_ _α_ _[M]_ _α_ _[M]_ _α_ [)] _[,]_ _σ_ _α_ _[M]_ _[∼]_ [Exponential(1)] _[.]_ 307 Similarly, for the trial random effects _β_ _j_ _[M]_ we assume _j_ _[M]_ _∼_ _N_ (0 _, σ_ _β_ _[M]_ _β_ _[M]_ _β_ [)] _[,]_ _σ_ _β_ _[M]_ _[∼]_ [Exponential(1)] _[.]_ 308 Also, these priors are weakly informative: the exponential distribution Exponential (1) has mean equal to 309 1. As a result, conditional on the values of _σ_ _α_ _[M]_ [and] _[σ]_ _β_ _[M]_ _[,]_ [the] [parameters] _[α]_ _i_ _[M]_ and _β_ _j_ _[M]_ are assigned a broad 310 prior, which nevertheless has a regularizing effect. 311 Finally, for the interaction random effect terms we assume _β_ _[M]_ _[,]_ [the] [parameters] _[α]_ _i_ _[M]_ _i_ _[M]_ and _β_ _[M]_ _[M]_ _ij_ _[∼]_ _[N]_ [(0] _[, σ]_ ( _[M]_ _α_ ( _α · β_ ) _[M]_ ( _α·β_ ) [)] _[,]_ _σ_ ( _[M]_ _α·β_ ) _[∼]_ [Exponential(1)] _[.]_ Priors used for parameters _t_ _[mid]_ 312 Priors used for parameters _t_ _ij_ and _r_ _ij_ can be found in the Supplementary Material accompanying 313 the article. Like for _M_ _[max]_ the is a ANOVA structure 313 the article. Like for _M_ _ij_ _[max]_, the key assumption is a two-way ANOVA structure _r_ _ij_ = exp( _θ_ _ij_ _[r]_ _ij_ _[r]_ [) = exp(] _[µ]_ _[r]_ [ +] _[ α]_ _i_ _[r]_ _j_ _[r]_ [+ (] _[α][ ·][ β]_ [)] _[r]_ _i_ _[r]_ [+] _[ β]_ _[r]_ _[r]_ _ij_ [)] _[,]_ (7) _t_ _[mid]_ _[mid]_ _ij_ = exp( _θ_ _ij_ _[mid]_ _ij_ _[mid]_ ) = exp( _µ_ _[mid]_ + _α_ _i_ _[mid]_ _i_ _[mid]_ + _β_ _[mid]_ _j_ _[mid]_ + ( _α · β_ ) _[mid]_ _ij_ _ij_ ) _._ (8) 314 Interpretation of the terms in these formulae is analogous to the one in (6). 315 For the measurement noise level we use the Half-Normal prior _σ_ _∼_ HalfNormal(0 _,_ 2) _._ 316 As the measurement errors _ε_ _ijk_ ’s are assumed to be additive and independent Gaussian, this leads to the 317 observation equation _y_ _ijk_ = _M_ ( _t_ _k_ ; _θ_ _ij_ ) + _ε_ _ijk_ = _M_ _ijk_ + _ε_ _ijk_ _,_ _ε_ _ijk_ _∼_ _N_ (0 _, σ_ ) _._ (9) 8 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.05.764527; this version posted November 6, 2026. The copyright holder for this](https://doi.org/10.64898/2026.11.05.764527) preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 318 Thus for a genotype _i_ and a trial _j_, the _y_ _ijk_ ’s are the longitudinal growth data. The model in (6), (7), (8) 319 and (9) is a GLMM with a non-standard link function for each of the DE parameters. The non-standard 320 link follows from a composition of the exponential function and the ODE solution (log-logistic function). 321 Normality of measurement errors is a reasonable assumption as a first approximation in the present setting. 322 In fact, daily APSIM wheat data were generated without noise, and therefore the unexplained residual 323 variation can be expected to be relatively small. The error standard deviation _σ_ quantifies the lack of fit 324 of the present model. The wide spacing of the data in time is expected to mitigate correlations between 325 residuals. See also Section 2.2 in Seber and Wild (1989). 326 Priors in Bayesian analyses can be set based on the experience of the user, knowledge of the problem at 327 hand, convenience, and regularization or numerical stability considerations. For data from plant breeding 328 experiments, the variation between trials is typically much larger than that between genotypes, while 329 _G × E_ is again smaller than variation between genotypes. We did not implement that knowledge in 330 the current analyses and used only weakly informative priors for our DE parameters, because we knew 331 that the data would be highly informative, while the prior effects would be limited. Nowadays, many 332 Bayesians are comfortable with using data-dependent priors; see, e.g., Section 9.5 in Gelman et al. (2021) 333 for the rationale behind default priors in the popular **R** package **rstanarm**, which is similar to the 334 **brms** procedure we used for our analyses in this paper, and cf. Gelman et al. (2017). The abundance of 335 data (specifically, the data for genotypes outside our subset of 48 genotypes) in our case offered some 336 possibilities for fine-tuning and experimentation with priors: for instance, we inspected where the observed 337 growth curves level off and centered the priors so that this range receives substantial prior mass. As a 338 qualification, we remark that we aimed at investigating dynamic phenotypes at medium-time resolution, 339 i.e., biweekly, as can occur in experiments with regular drone flights. The number of observations per 340 growing season for such experiments is maximally around 10 to 15. This number restricts the complexity 341 of the DE system that can be successfully fitted to such data. 342 **2.6** **Parameter** **identifiability** Although in Section 2.5 the logarithms of the parameters of the differential equation, i.e. _θ_ _[M]_ _ij_ _[r]_ [and] _[ θ]_ _ij_ _[tmid]_ _ij_ _[, θ]_ _[r]_ 343 _ij_ _ij_ _ij_ _,_ 344 are identifiable, their additive decompositions into the sums of the main genotypic and environmental 345 effects and their interactions are not unique. This is well-known for ANOVA-type models. Additional 346 restrictions must be imposed to talk about genotypic, environmental, and _G × E_ effects. Thus, once 347 the MCMC samples of the parameters have been collected, in the post-processing step we applied the 348 sum-to-zero restrictions on _α_ _i_ _[p]_ _[,]_ _[β]_ _j_ _[p]_ [and] [(] _[α][ ·][ β]_ [)] _[p]_ _ij_ _[,]_ [with] _[p]_ [being] [one] [of] _[M]_ [,] _[t]_ _[mid]_ [and] _[r.]_ [See,] [e.g.,] [Kruschke] 349 (2014), Section 20.1.1, Gelman and Hill (2006), Section 19.4, and Ogle and Barber (2020) for details. Main 350 effects and interactions refer to parameters in this new parametrization, and are unambiguously defined. 351 In principle, the new parameters need to receive other names, such as _α_ ˜︁ _i_ _[p]_ _[,]_ _[β]_ [˜︁] _j_ _[p]_ [and] [(] _α_ [ ˜︃] _· β_ ) _[p]_ _ij_ _[.]_ [However,] [for] _i_ _[p]_ _[,]_ _[β]_ _[p]_ _j_ _[p]_ [and] [(] _[α][ ·][ β]_ [)] _[p]_ _j_ _[p]_ [and] [(] _α_ [ ˜︃] _· β_ ) _[p]_ _i_ _[p]_ _[,]_ _[β]_ [˜︁] _j_ _[p]_ 351 In principle, the new parameters need to receive other names, such as _α_ ˜︁ _i_ _[p]_ _[,]_ _[β]_ _j_ _[p]_ [and] [(] _α · β_ ) _[p]_ _ij_ _[.]_ [However,] [for] 352 notational simplicity, in the rest of the paper we will keep the old notation _α_ _i_ _[p]_ _[,]_ _[β]_ _j_ _[p]_ [and] [(] _[α][ ·][ β]_ [)] _[p]_ _ij_ _[,]_ [but] [we] 353 will mean the parameters in the new parameterization. We also remark that priors in Bayesian mixed 354 models help with weak identifiability of parameters. In fact, the conclusions in this paper would have 355 stayed largely unchanged, had we used in our analysis the raw random effects instead of the adjusted ones. _i_ _[p]_ _[,]_ _[β]_ _[p]_ _j_ _[p]_ [and] [(] _[α][ ·][ β]_ [)] _[p]_ 356 **2.7** **Interpretation** **of** **genotypic,** **environmental** **and** _G × E_ **parameters** 357 Interpretations of dynamic behavior shown by individual genotypes and induced by individual envi358 ronments can be based on estimates for individual effects and their dispersion parameters building on 359 the interpretation of similar parameters in traditional quantitative genetic models for continuous traits. 360 To be more specific, we will build such interpretations on the posterior distributions for genotypic and 361 environmental main effects and _G × E_ effects, i.e., _α_ _i_ _[p]_ [,] _[β]_ _j_ _[p]_ [and] [(] _[α][ ·][ β]_ [)] _[p]_ _ij_ [,] [respectively,] [with] _[p]_ [being] [one] [of] 362 _M_, _t_ _[mid]_ and _r_, together with the corresponding dispersion parameters _σ_ _α_ _[p]_ [,] _[σ]_ _β_ _[p]_ [,] [and] _[σ]_ ( _[p]_ _α.β_ ) [.] [Inspection] [of] 363 the posterior distributions of the dispersion parameters will inform us about the relative importance of 364 genotypic and environmental main effects as well as that of _G × E_ . To appreciate the relative magnitudes 365 of _G × E_ effects versus genotypic main effects for the dynamic parameters, we will look at the posteriors 366 for the ratios of the _G × E_ standard deviations to the standard deviations for genotypic main effects. 367 For dynamic parameters with sizeable dispersion for _G × E_, individual interaction effects, ( _α · β_ ) _[p]_ _ij_ [,] 368 can be studied via their credible intervals to verify whether they contain zero. Questions concerning 369 genotypic and environmental characterizations can be addressed by counts on non-null _G × E_ effects. 370 Genotypes can be characterized on the percentage of times they show non-null _G × E_ effects across all 371 trials or subsets of trials. This percentage can be interpreted as a measure of stability. The four locations _i_ _[p]_ [,] _[β]_ _[p]_ _j_ _[p]_ [and] [(] _[α][ ·][ β]_ [)] _[p]_ _α_ _[p]_ [,] _[σ]_ _[p]_ _[p]_ [and] _[σ]_ _[p]_ _β_ [,] ( 9 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.05.764527; this version posted November 6, 2026. The copyright holder for this](https://doi.org/10.64898/2026.11.05.764527) preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 372 that were included in the simulation data represent different growing regions within Australia. Therefore, 373 it may be of interest to see to which extent the incidence of non-null _G × E_ effects differs between those 374 regions. A low incidence may point to a non-discriminative environment, whereas a high incidence would 375 point to a discriminative environment. 376 Posterior distributions of environmental and/or _G × E_ effects can be used for environmental character377 ization to create so-called environmental scenarios or envirotypes (Xu, 2016; Resende et al., 2021; Chen 378 et al., 2026; Fernandez et al., 2026). Typically in environmental characterization, trials are clustered into 379 groups, with the groups representing a distinguishable type of growing environment. Envirotypes can 380 serve to develop more efficient breeding strategies. We illustrate how to use posterior medians for the 381 environmental main effect parameters _β_ _j_ _[M]_ [,] _[β]_ _j_ _[mid]_, and _β_ _j_ _[r]_ [as] [inputs] [for] [an] [environmental] [characterization.] 382 These posterior medians were standardized and collected in a vector of length three. These vectors were 383 clustered into three groups using a _K_ -means algorithm (MacQueen, 1967). Choosing three or four clusters 384 is supported by relevant literature on envirotypes (Fernandez et al., 2026). The results of the _K_ -means 385 clustering depend on the random initialization of the algorithm. With three groups, the inferred clusters 386 were stable over multiple runs of the algorithm. 387 For plotting purposes, we used for some figures medoid representations (Van der Laan et al., 2003). For 388 environments, this resulted in real environments that were at the center of the space spanned by estimates 389 for the environmental main effects in the GLMMs for the three parameters. The medoid environment can 390 then be thought of as a typical environment within a class of environments. The same technique was also 391 used for plotting genotypic information. _j_ _[M]_ [,] _[β]_ _j_ _[mid]_ _j_, and _β_ _j_ _[r]_ 392 **2.8** **Software** **implementation** 393 We implemented our Bayesian approach in **Stan**, specifically using the **R** package **brms** . **Stan** is a 394 state-of-the-art platform for statistical modeling and high-performance statistical computation. The 395 principal engine to perform Bayesian inference in **Stan** is the No-U-Turn-Sampler (NUTS), a cutting-edge 396 variant of the Hamiltonian Monte Carlo (Hoffman and Gelman, 2014; Betancourt, 2017; Betancourt 397 et al., 2017). The **R** package **brms** uses **Stan** as its backend, but employs the formula syntax similar 398 to the popular **lme4** package. This results in a familiar and simple interface for performing regression 399 type analyses. See B¨urkner (2017, 2018, 2021) for additional information. For a general overview and 400 comparison of the modern software for Bayesian inference, see Strumbelj [ˇ] et al. (2024). 401 **brms** automatically evaluates the popular _R_ [ˆ] convergence diagnostic for Markov chains it generates 402 (Gelman and Rubin, 1992; Vehtari et al., 2021), and issues a warning when its values are above the 403 commonly used threshold 1 _._ 05 _._ See also Table 5 in the Supplementary Material for values of _R_ [ˆ], which are 404 all below the threshold 1.05. In order to verify convergence, we also visually examined the trace plots 405 and autocorrelation functions for our Markov chains to the extent that this is possible given the large 406 number of parameters in the model (Gabry et al., 2025). The plots for an arbitrary selection of the model 407 parameters are presented in the Supplementary Material in Figures 11 and 12. From Table 5 we also see 408 that the lowest effective sample sizes occur for the global intercept parameters, which are weakly identified 409 from the data and are not of inferential interest. On the other hand, the standard deviation parameters 410 have bulk and tail effective sample sizes above 400 _,_ which allows for accurate simulation-based inference 411 in the present context. 412 To handle the **brms** output, we mainly employed the **R** package **tidybayes** (Kay, 2022), and other 413 capabilities offered by **R** (R Core Team, 2022). For plotting in **R** we relied on **ggplot2** (Wickham, 2016), 414 as well as several other packages building upon and extending **ggplot2** . 415 In our implementation in **brms**, we used an explicit form of the log-logistic function given in (4) . 416 For future extensions, we note that availability of an analytic solution to a differential equation is not 417 a requirement for our Bayesian approach: it also works when the analytic solution is replaced with a 418 numerical solution. In that case, instead of **brms**, the model may be implemented directly in **Stan** ; cf. 419 Chatzilena et al. (2019) and Grinsztajn et al. (2021). Alternatively, additional **Stan** code can be injected 420 into the **brms** model specification. Fitting, however, is faster with an analytic solution. As for the 421 effect of the approximation error on posterior inference resulting from numerical integration, the impact 422 stays bounded, and is negligible when the step size of the ODE solver is chosen to be small enough; see 423 Donnet and Samson (2007) for a mathematically precise statement. **Stan** implements a number of reliable 424 adaptive ODE solvers (Stan Development Team, 2023, Chapter 13). For early contributions to Bayesian 425 fitting of hierarchical differential equation models, see Bois et al. (1996); Gelman et al. (1996); Wakefield 426 (1996). These early contributions rely on Gibbs sampling (Geman and Geman, 1984; Gelfand and Smith, 427 1990). However, the Gibbs sampler is known to suffer from poor mixing in statistical problems with 10 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.05.764527; this version posted November 6, 2026. The copyright holder for this](https://doi.org/10.64898/2026.11.05.764527) preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 428 highly correlated parameters in the posterior (De Jong and Shephard, 1995; Lindsten and Sch¨on, 2013). 429 This is inherent in the sampler because of its sequential updating of single individual parameters in any 430 given iteration. This slow convergence hampers the scalability of the Gibbs sampler to high-dimensional 431 models with complex posterior geometry, such as those defined by differential equations (Girolami, 2008; 432 Calderhead and Girolami, 2009), and especially their hierarchical extensions. ## 433 3 Results 434 **3.1** **Goodness** **of** **fit** **and** **characterization** **of** **sources** **of** **variation** 435 To assess how well our model fitted the ‘observed’ longitudinal data for individual genotypes in individual 436 environments, we looked at plots of observed (Figure 2) and fitted (Figure 3) data. The observations were 437 created by sampling at medium-time resolution biomass data from APSIM wheat simulations. There were 438 four locations with 31 years of trials, in total 124 environments, for a set of 48 genotypes. In Figure 2, for 439 each location we present the so-called medoid trial, a trial that can be thought of as best representing the 440 group (location) to which it belongs. For such a medoid trial, we show the daily data for all genotypes. 441 The colors of the genotypes are consistent in Figures 2 and 3. Moreover, to give an idea of the amount 442 of uncertainty in the obtained fits, in Figure 4 we show 95% Bayesian pointwise prediction bands for 443 the medoid genotype in the medoid year (trial) at each of the 4 locations. As can be readily seen, most 444 observations fall within the prediction bands. 445 The log-logistic fits were good in that they follow the observed data closely for all genotypes and trials, 446 and for most predictive purposes in this setup the log-logistic fits are probably good enough. Still, the 447 observed data contain small bumps, deviations from the smooth log-logistic fits, which are very likely 448 caused by short-term meteorological events, such as elevated temperatures and drought. Such bumps 449 cannot be incorporated into the log-logistic model, but can show up in residual plots, see Figure 13. 450 Figure 2 contains genotypic curves that hardly intersect, but there are clear differences between 451 locations, between years within locations (trials) and between genotypes. Narrabri and Yanco attained 452 higher biomasses than Merredin and Emerald, with Emerald having the lowest biomasses. The amount of 453 variation between years within locations seems comparable across the four locations, see Figure 5. 454 Inspection of the posterior standard deviations for environmental and genotypic main effects and 455 _G × E_ effects in Figure 6 shows that the dynamic parameters for asymptote, midpoint and slope follow 456 the standard partitioning of phenotypic variation with environmental main effects being responsible for 457 the majority of the phenotypic variation, followed by genotypic main effects and finally _G × E_ effects. 458 The greatest variation was found in the asymptote, followed by the midpoint, while the slope was clearly 459 less variable in this wheat panel. _G × E_ variation is best appreciated in relation to genotypic main effect 460 variation. Figure 7 presents ratios of _G × E_ standard deviation to genotype standard deviation. For the 461 asymptote, this ratio was between 0 _._ 6 and 0 _._ 9 _,_ while for the midpoint it was between 0 _._ 25 and 0 _._ 35 _._ For 462 the slope, the same ratio had a value between 0 _._ 0 and 0 _._ 1 _._ The _G × E_ appeared clearly for the asymptote, 463 but for the midpoint and slope parameter it had little to no influence. Therefore, the rough biomass 464 growing pattern for wheat seems predictable from additive GLMMs for at least midpoint and slope. That 465 is, no or hardly any _G × E_ plays a role for midpoint and slope parameters in the biomass growth functions. 466 For the asymptote, the maximum biomass, _G × E_ does play a role, and further study is required to find 467 out how to arrive at a predictive model for biomass taking into account the _G × E_ for maximum biomass. 468 **3.2** _G × E_ **interactions** 469 Interactions can be further studied via credible intervals for interaction terms ( _α · β_ ) _[p]_ _ij_ [,] [and] [verifying] 470 whether zero is contained within these intervals or not. Table 2 summarizes the results and confirms the 471 impression that with the log-logistic model, most of the interaction is to be discovered for the asymptote, 472 _M_ _ij_ _[max]_ . For this DE parameter, _G × E_ did occur in about a third of all _G × E_ combinations across the full 473 set of trials. More detail is provided by the number or percentage of _G × E_ interactions per location in 474 Table 3. For parameters _M_ _ij_ _[max]_, there is some variation between locations in the frequencies of interactions, 475 with _G × E_ being more frequent in Yanco and less frequent in Merredin. Further zooming on _G × E_ 476 interactions, we inspect the non-zero interactions in the trials over time (years) for each of the locations. 477 In Figure 8, we again see that Yanco in general induces more _G × E_ than the other three locations, with 478 Merredin inducing the lowest frequency of _G × E_ . We also see that there is a lot of variation in the 479 frequency of non-null _G × E_ between years within locations. 11 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.05.764527; this version posted November 6, 2026. The copyright holder for this](https://doi.org/10.64898/2026.11.05.764527) preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. _[M]_ _ij_ [,] [(] _[α][ ·][ β]_ [)] _[mid]_ _ij_ Table 2: Frequencies of credible intervals for interaction terms ( _α · β_ ) _[M]_ _ij_ and ( _α · β_ ) _[r]_ Table 2: Frequencies of credible intervals for interaction terms ( _α · β_ ) _ij_ [,] [(] _[α][ ·][ β]_ [)] _ij_ and ( _α · β_ ) _[r]_ _ij_ [that] exclude zero, i.e., that do indicate _G × E_ . _[M]_ _ij_ **(** _α · β_ ) _[mid]_ _ij_ **(** _α · β_ ) _[M]_ _ij_ _ij_ **(** _α · β_ ) _[r]_ _ij_ _ij_ Credible interval excludes zero 2173 91 0 Percentage of intervals without zero 36 _._ 5% 1 _._ 5% 0% _[M]_ _ij_ [,] [(] _[α][ ·][ β]_ [)] _[mid]_ _ij_ Table 3: Frequencies of credible intervals for interaction terms ( _α · β_ ) _[M]_ _ij_ and ( _α · β_ ) _[r]_ Table 3: Frequencies of credible intervals for interaction terms ( _α · β_ ) _ij_ [,] [(] _[α][ ·][ β]_ [)] _ij_ and ( _α · β_ ) _[r]_ _ij_ [that] exclude zero, i.e., that do indicate _G × E_ . Frequencies are grouped by location. _[M]_ _ij_ **(** _α · β_ ) _[mid]_ _ij_ _ij_ **(** _α · β_ ) _[M]_ _ij_ _ij_ **(** _α · β_ ) _[r]_ _ij_ Emerald 534 0 0 Merredin 422 15 0 Narrabri 586 41 0 Yanco 631 35 0 Total 2173 91 0 480 When we apply the same logic to the genotypes, we can compare genotypes for the number of non-null 481 _G × E_ effects across trials. We can interpret this number as an indicator for the genotypic stability of DE 482 parameters. The lower the frequency or number of non-null interactions, the higher the stability of the 483 genotype. Figure 9 offers a view on substantial genotypic variation for non-zero interactions for _G × E_ in 484 the asymptote, ( _α · β_ ) _[M]_ _ij_ [.] 485 **3.3** **Environmental** **characterization** 486 The parameters in the GLMMs for the DE parameters can be used to characterize genotypes and 487 environments. For the environments, the environmental main effect parameters _β_ _j_ _[M]_ [,] _[β]_ _j_ _[mid]_, and _β_ _j_ _[r]_ [should] 488 reflect the local growing conditions in ways comparable to physical and physiological measurements. Upon 489 collecting these parameters in three-dimensional vectors, we used the _K_ -means clustering algorithm on 490 these environmental descriptors. The _K_ -means solution was stable for solutions with three clusters over 491 multiple runs with different initializations. On inspection of these clusters, it turned out that Emerald 492 and Merredin trials each dominated one cluster, while a third cluster was a combination of especially 493 Narrabri and to a lesser extent Yanco trials. Envirotypes facilitate the efficient modeling of _G × E_ by 494 partitioning the initial _G × E_ in a part explained by _G × Envirotype_ interaction and a residual part. For 495 the APSIM wheat data here, because the envirotypes were dominated by differences between locations, a 496 description of _G × E_ for the asymptote by introducing a _G × Envirotype_ term was about as effective as 497 introducing a _G × Location_ term, and a comparable residual _G × E_ remained. _j_ _[M]_ [,] _[β]_ _j_ _[mid]_ _j_, and _β_ _j_ _[r]_ 498 **3.4** **Genotypic** **characterization** 499 In Bustos-Korts et al. (2019b), genotype-specific values for a panel of 199 genotypes were generated for 12 500 APSIM parameters that determine phenology, capture of environmental resources, resource use efficiency, 501 and biomass partitioning. APSIM parameters were created from a genetic basis consisting of 300 QTLs 502 whose effects were sampled from gamma distributions with parameters following estimates obtained from 503 analysing real wheat data. 504 We wanted to characterize the genotypes in our hierarchical analysis by the estimated genotypic 505 main effects in the GLMMs for the DE parameters. To this end, we compared the estimates for these 506 genotypic effects with APSIM parameters used to generate the data by Bustos-Korts et al. (2019b). 507 APSIM parameters that changed over time were not taken into account. The selected APSIM parameters 508 are summarized in Table 4, which is constructed as a subset of Table 2 in Bustos-Korts et al. (2019b). 509 We computed correlations between each of the DE parameter _α_ _i_ _[M]_ _[, α]_ _i_ _[mid]_ and _α_ _i_ _[r]_ [(estimated by posterior] 510 medians) on the one hand, and the individual APSIM parameters, on the other. We also obtained the 511 asymptotic 95% confidence intervals based on the Fisher z-transformation using the standard **R** command 512 **cor.test** . The results are presented graphically in Figure 10. 513 Several APSIM parameters show significant correlations with the inferred parameters of our Bayesian 514 hierarchical model. The maximum biomass seems to be strongly negatively related to the lower limit 515 for water uptake ( `ll` ~~`m`~~ `odifier` ), whereas that relationship is positive for the slope parameter. Higher _i_ _[M]_ _[, α]_ _i_ _[mid]_ _i_ and _α_ _i_ _[r]_ 12 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.05.764527; this version posted November 6, 2026. The copyright holder for this](https://doi.org/10.64898/2026.11.05.764527) preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. Table 4: Genotype-specific parameters from APSIM that were used in this study for comparison purposes. **Parameter** **APSIM** **name** **Unit** Fr. of leaves senescing main stem node `fr` ~~`l`~~ `f` ~~`s`~~ `en` ~~`r`~~ `ate` [ ] Grains/gram of stem at flowering `grains` ~~`p`~~ `er` ~~`g`~~ `ram` ~~`s`~~ `tem` grains _g_ _stem_ _[−]_ [1] Lower limit for water uptake `ll` ~~`m`~~ `odifier` [ ] Maximum grain size `max` `grain` ~~`s`~~ `ize` mg Sensitivity to photoperiod `photop` ~~`s`~~ `ens` [ ] Potential grain filling rate `potential` ~~`g`~~ `rain` ~~`f`~~ `illing` ~~`r`~~ `ate` g grain _[−]_ [1] day _[−]_ [1] Thermal time to reach floral initiation `tt` ~~`f`~~ `loral` ~~`i`~~ `nitiation` _◦_ C d Sensitivity to vernalization `vern` `sens` [ ] 516 slopes were associated with larger thermal times till floral initiation ( `tt` ~~`f`~~ `loral` ~~`i`~~ `nitiation` ) and stronger 517 photoperiod sensitivity ( `photop` ~~`s`~~ `ens` ). The midpoint time did not show any obvious correlation with 518 APSIM parameters. ## 519 4 Discussion 520 **4.1** **Incorporating** **phenotyping** **information** **into** **mixed** **models** **for** **genomic** 521 **prediction** 522 In recent years, automatic monitoring of longitudinal phenotyping data has become an integral part of 523 many plant breeding programs and experiments. Incorporating this kind of information into breeding 524 strategies and decisions is a topic of current research. From a prediction perspective, genomic prediction 525 models can be improved by combining genomic and phenotyping data into genotype-to-phenotype models. 526 An important objective of such integrated models is to increase the predictive ability of within-trial 527 phenotypes at later time points from genomic information and phenotyping information at earlier time 528 points. In other words, the target phenotypes, such as yield, phenological traits, and yield components, 529 are forecasted from early season information. The chosen modeling framework is often that of multi-trait 530 mixed models and random regressions (Baba et al., 2020; Campbell et al., 2018; Momen et al., 2019; Sun 531 et al., 2017, 2019). These models can contain complex formulations for the variance-covariance structures 532 dealing with the correlations between traits over time and the incorporation of one or more kernels on the 533 genotypes to represent genomic relations, but they do not include explicit model terms or structures that 534 represent physiological mechanisms or insights. Furthermore, the above papers leave open how genomic 535 and phenotyping information should be combined to predict yield in future environments under _G × E_ . 536 Jarqu´ın et al. (2014) presented a mixed model framework for genomic prediction between trials in the 537 presence of _G_ _×_ _E_ based on the definition of kernels (relationship matrices) for genotypes and environments. 538 A generalization of this approach addresses _G × E_ while integrating genomic and phenotyping information 539 on genotypes with enviromic information on environments (Adak et al., 2023; Robert et al., 2022). 540 Although the emphasis for this new type of integrated mixed models is on between-trial prediction in 541 place of within-trial prediction, the representation of the additional phenotyping information inside the 542 prediction models ignores most of the physiological principles that make yield the outcome of growth 543 and developmental processes. In this approach, the dynamic environmental inputs that are typical for 544 physiological models are not used. The totality of the physiological information is reduced to a relationship 545 matrix between genotypes. **–** 546 **4.2** **Crop** **growth** **model** **whole** **genome** **prediction** 547 In stark contrast to the mixed models above that incorporated phenotyping information without reference 548 to physiological models, an alternative approach, called crop growth model - whole genome prediction, 549 builds on elaborate physiological models for prediction in the expectation that such an approach will lead 550 to better between-trial predictions of yield under _G × E_ (Technow et al., 2015; Messina et al., 2018). The 551 CGM-WGP approach predicts the end-point trait yield together with dynamical yield components for 552 new genotypes and environments using marker profiles and dynamic environmental data as inputs. The 553 marker profiles serve to provide a molecular basis to a limited number of physiological traits. Marker 554 effects are estimated following a Bayesian procedure. The CGM-WGP represents a powerful, general and 555 attractive predictive framework founded on prior physiological knowledge. In its initial formulations, the 13 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.05.764527; this version posted November 6, 2026. The copyright holder for this](https://doi.org/10.64898/2026.11.05.764527) preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 556 CGM-WGP approach does not use dynamic observations on plant growth during the season and has to 557 rely fully on the correctness of the pre-specified parameters and relations inside the crop growth model 558 that converts marker profiles and environmental inputs into yield outputs. 559 Crop growth models tend to have many parameters and it is impossible to estimate all of the parameters 560 in crop growth models at a genotype-specific level. When using a CGM-WGP approach, a decision has to 561 be made about which parameters to allow to vary with genotype and which to fix for the crop species as a 562 whole. In the CGM-WGP papers by Technow et al. (2015), Cooper et al. (2016), Messina et al. (2018) 563 and Jighly et al. (2023) only a few yield components are genotype specific and a function of marker inputs. 564 Although the CGM-WGP approach is dynamic for the inputs from the environment and the outputs of 565 the yield components, it does not use dynamic inputs from phenotyping and phenomic devices. Such 566 inputs may contain important information on the occurrence of biotic and abiotic stresses. 567 **4.3** **Our** **hierarchical** **dynamic** **model** **proposal** **to** **function-valued** **modeling** **of** 568 **longitudinal** **traits** **in** **multiple** **trials** 569 The Bayesian framework we developed in this paper provides a flexible framework to model a longitudinal 570 trait for a set of genotypes across multiple environments, taking into account _G × E_ . Our approach is 571 unique because it fits the growth curves simultaneously for all included genotypes and environments. 572 Therefore, propagation of errors and tests about differences between growth trajectories can be performed 573 within a well-defined context, namely that of the Bayesian framework that relies for inference on the 574 properties of the posterior distribution. The same framework contains the tools for constructing prediction 575 intervals along the fitted growth trajectories. Each trajectory is determined by a limited number of 576 biologically interpretable DE parameters whose expectation is itself a GLMM with main effects for 577 genotype and environment, as well as _G × E_ . 578 In our approach, we followed a function-valued trait modeling perspective to growth trajectories as 579 propagated by Wang et al. (2014), Gomulkiewicz et al. (2018) and Baker and Wang (2021). Within that 580 framework, we interpret the DE parameters as phenotypic traits themselves, as traits underlying the 581 dynamics of the yield component or secondary trait that we are studying. Analogously to a quantification 582 of the sources of variation in the context of linear mixed models in quantitative genetics, the contributions 583 of genotype, environment, and _G × E_ to the total ’phenotypic variation’ for the DE parameter can be 584 assessed within the Bayesian framework by inspecting the posterior distributions for the corresponding 585 dispersion parameters. Their magnitudes reflect their relative importance. For classical phenotypes such 586 as yield, biomass and phenology, the environment is by far the most important source of variation, usually 587 followed by the genotype when the environmental conditions are not too wildly different between trials 588 and finally the _G × E_ . Exactly that order of importance also occurred for the DE parameters in our 589 model for biomass, for all three the DE parameters of asymptote, midpoint, and slope. Comparison of 590 dispersion parameters between the three DE parameters identified the asymptote as the parameter with 591 the strongest _G × E_, that is, with the highest ratio of _G × E_ to genotype variation. For the midpoint and 592 slope, _G × E_ was small to absent, which means that the dynamics of biomass followed an additive model 593 for these two parameters. 594 **4.4** **Alternative** **function-valued** **modeling** **approaches** **for** **longitudinal** **traits** 595 A function-valued approach to plant growth for populations of genotypes was also proposed by Baker et al. 596 (2018a,b). For leaf width and length in _Brassica_ _rapa_, logistic models were fitted on growing degree days 597 per plant. The growth parameters were then extracted and subsequently further analyzed in a Bayesian 598 hierarchical analysis of variance (ANOVA) model containing a population, genotypic and plant level. 599 Both Baker et al. (2018a) and Baker et al. (2018b) tested plant growth under two plant densities that 600 were interpreted as different environmental conditions. A test for genotype by density interactions was 601 performed and on failure to reject the null hypothesis of no interaction, it was concluded that the growth 602 of a genotype as modeled under one density could be used to predict the growth of the same genotype 603 under an alternative density. In our approach, in contrast to Baker et al. (2018a) and Baker et al. (2018b), 604 a single inferential framework was used for both modeling the dynamics of plant growth and expressing 605 the dynamic parameters in terms of genotype, environment, and _G × E_ . A test of _G × E_ for individual DE 606 parameters followed straightforwardly from inspection of the posterior for the corresponding dispersion 607 parameter. 608 Baker et al. (2018a) extended their Bayesian ANOVA model for the logistic parameters by making the 609 genotypic means for these parameters dependent on photosynthetic capacity. In this way, they wanted 14 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.05.764527; this version posted November 6, 2026. The copyright holder for this](https://doi.org/10.64898/2026.11.05.764527) preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 610 to test whether carbon limitations were responsible for genotypic differences as a factor affecting the 611 growth dynamics of individual genotypes. The photosynthetic capacity can be measured in different 612 ways, some of them being very time consuming others delivering approximations via remote sensing. The 613 insertion of approximate physiological parameters from phenotyping devices or phenomic measurements 614 to describe variation in growth parameters offers an attractive and powerful modeling perspective. Similar 615 tests for DE parameters being a function of certain genotype-dependent limitations or capacities can 616 naturally be implemented in our framework by introducing genotypic covariates into the GLMMs for the 617 DE parameters and inspecting the posterior distribution for the slope with respect to such a covariable. 618 Baker et al. (2018b) identified quantitative trait loci (QTLs) for dynamic parameters, after estimation 619 of those dynamic parameters in logistic growth models. They follow up by creating QTL-based predictions 620 for the full growth trajectories. Again, this is a multi-step approach. Within our Bayesian hierarchical 621 model for plant growth, we can in principle directly predict the dynamical response from marker profiles 622 by replacing the genotypic main effects in the GLMMs for the growth parameters by penalized regressions 623 on marker profiles. With a double exponential prior on the marker effects this would lead to QTL mapping 624 via LASSO (Least Absolute Shrinkage and Selection Operator) regression (Li and Sillanp¨a¨a, 2012), while 625 with a normal prior we would arrive at genomic prediction using a ridge regression. 626 In the field of animal breeding, Varona et al. (1999) use the Gibbs sampler (or, rather, the Metropolis627 within-Gibbs sampler) with an explicit form of the growth curve, specifically the Von Bertalanffy growth 628 curve, within the mixed effects modelling framework. Their approach is not generalisable to the differential 629 equation models that do not admit explicit solution. The Gibbs sampler also does not scale well in such 630 complex nonlinear models. In fact, Onogi et al. (2019), who use a similar computational approach as 631 Varona et al. (1999), are constrained to run their Gibbs sampler for millions of iterations (Onogi et al., 632 2019, p. 5). This is indicative of slow mixing of the sampler. In contrast, in order to efficiently explore the 633 posterior distribution, a few thousand iterations typically suffice with the NUTS sampler (Hoffman and 634 Gelman, 2014; Betancourt, 2017; Betancourt et al., 2017), which is the computational algorithm of **brms** . 635 The mathematical model in Onogi (2020) is similar to that in Varona et al. (1999); Onogi et al. (2019) 636 and does not envision a governing differential equation without a closed-form solution. The computational 637 approach they use is a hybrid method from Onogi et al. (2015). This combines Gibbs sampling with 638 coordinate ascent variational inference (CAVI). However, even though in general faster than MCMC 639 sampling, CAVI is closely related to Gibbs sampling and suffers from the same general drawbacks (Blei 640 et al., 2017). A computational bottleneck of the Bayesian approach to inference in DE models lies in the 641 need for their repeated numerical integration in likelihood evaluations. Already for a medium complexity 642 model, this is not realistically feasible when the computational method requires millions of iterations 643 until convergence (Voit and Almeida, 2004, Appendix 1). Note also that while the MCMC methods are 644 guaranteed (under idealised conditions) to recover the posterior distribution as the number of iterations 645 grows indefinitely (Roberts and Rosenthal, 2004), this is not the case with variational inference, which 646 involves a non-vanishing approximation error even in the limit. 647 **4.5** **Outlook** 648 The current version of our hierarchical Bayesian model for longitudinal traits in multiple trials is promising 649 as a starting point for a framework for the prediction of yield under _G × E_ . However, a number of issues 650 require attention. First of all, from the candidate DE models described in Van Voorn et al. (2023) we 651 chose for this paper the simplest one, the logistic model, which was slightly modified, for computational 652 reasons, into a log-logistic model. The log-logistic model does not contain an environmental driver, as in 653 the irradiance model or the temperature model of Van Voorn et al. (2023). Neither does the log-logistic 654 model contain an environmental state variable, as in the water model of Van Voorn et al. (2023). The 655 log-logistic model has an explicit solution that turns our hierarchical dynamic model into a hierarchical 656 regression model. Therefore, estimation procedures are obviously easier to perform. Nevertheless, the 657 strong point of our approach remains that parameters for growth processes can explicitly be modeled as 658 functions of genotypic, environment and _G×E_ terms. The usefulness of parameters of yield components for 659 incorporation into predictive models for yield will firstly depend on the extent to which these parameters 660 are free of _G × E_ . Even when parameters for yield components do show _G × E_, they can still be useful for 661 yield prediction as long as the _G × E_ can be modeled in terms of environmental covariables, that is, the 662 _G × E_ term for the rate parameter is written as the product of a genotypic sensitivity times a predictable 663 environmental condition. The situation here is comparable to the use of factorial regression in models for 664 _G × E_ in target traits (Van Eeuwijk et al., 2016; Avagyan et al., 2025). 665 When our DE system does not have an explicit analytical solution, our hierarchical Bayesian framework 15 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.05.764527; this version posted November 6, 2026. The copyright holder for this](https://doi.org/10.64898/2026.11.05.764527) preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 666 will have to solve the DE system for each combination of genotype and environment numerically together 667 with estimating the GLMM parameters for the DE parameters. This is possible within our modeling 668 framework that relies on the Bayesian modeling language **Stan** and we will present results in follow-up 669 papers. Computing resources will need to be increased when there is no explicit solution for the DE 670 system. For the analyses presented in this paper, we kept the playing field small as the objective was to 671 investigate the principles of a hierarchical Bayesian approach to longitudinal modeling. A first next step 672 will consist of upscaling our modeling facilities to medium sized panels of genotypes (100s) and sets of 673 environments (25–100). 674 An ultimate aim of longitudinal traits modeling in phenotyping experiments is to arrive at a predictive 675 model that allows us to predict yield under _G × E_ . Different DE systems can be tried to see whether 676 formulations can be found in which parameters have no or limited _G × E_ or in which _G × E_ depends on 677 predictable environmental covariables. Further, the genotypic terms, intercepts, and slopes in the GLMMs 678 for the DE parameters should be predictable from marker profiles, like in standard genomic prediction 679 approaches. For the prediction of genotypic intercepts and slopes for DE parameters, an additional 680 layer can be added to the Bayesian hierarchy. An alternative is a two-step approach in which these 681 genotypic intercepts and slopes are predicted in a separate genomic prediction step. An interesting aid in 682 identifying candidate DE systems for modeling longitudinal traits is the SINDy-PI algorithm (Kaheman 683 et al., 2020) that tries to identify sparse representations of dynamic systems. An illustrative application 684 of this approach in the context of genotype-specific crop growth modeling is presented by Verburg et al. 685 (2026). We feel that our hierarchical Bayesian approach to longitudinal trait modeling opens up new 686 perspectives to predict yield under _G × E_ . This approach provides an additional tool for breeders and 687 geneticists in conjunction with the mixed kernel-based models following the strategy outlined by Jarqu´ın 688 et al. (2014) and the CGM-WGP strategy introduced by Technow et al. (2015) and Messina et al. (2018). ## 689 Conflict of interest statement 690 Authors ST, NF, RM and CM were employees of Corteva Agriscience when this manuscript was conceived. 691 The remaining authors declare that the research was conducted in the absence of any commercial or 692 financial relationships that could be construed as a potential conflict of interest. ## 693 Author contributions 694 Modelling **G** _×_ **E** _×_ **M** was conceived between Corteva Agriscience and Wageningen University & Research 695 – Biometris, with groups led by CM and FV, respectively. Conceptualization and development was done 696 by all authors. The APSIM wheat simulation data were created by DB. Model fitting was done by SG. 697 The first version manuscript writing and preparation was done by SG, GV and FV. Intermediate editing 698 of the manuscript was done by all authors. SG, GV, MB and FV did the final manuscript editing. All 699 authors contributed to the article and approved the submitted version. ## 700 Funding 701 This research was supported by Predictive Agriculture, Research & Development at Corteva Agriscience 702 (2019–2020). 703 FV worked on this paper partly in the context of the EU-INVITE project, which is part of the 704 Horizon 2020 Framework Programme of the European Union, under Grant Agreement 817970 ( `[https:](https://www.h2020-invite.eu/)` 705 `[//www.h2020-invite.eu/](https://www.h2020-invite.eu/)` ). 706 FV, MB and SC contributed to this paper partly within the framework of the GRDC INVITA UQ-WUR 707 and GRDC AAGI UQ-WUR Collaborative Research Agreements. 708 SG, MB and GV acknowledge funding from the WUR Investment Theme on Data-Driven Discoveries 709 in a Changing Climate (D3-C2). 710 DB contributed to this paper thanks to the funding of Subvenci´on a la Instalaci´on en la Academia, 711 ANID 85220037. 16 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.05.764527; this version posted November 6, 2026. The copyright holder for this](https://doi.org/10.64898/2026.11.05.764527) preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. ## 712 Acknowledgments 713 Acknowledgment is made to the APSIM Initiative, which takes responsibility for quality assurance and a 714 structured innovation program for APSIM’s modelling software, which is provided free for research and 715 development use (see `[https://www.apsim.info/](https://www.apsim.info/)` for details). ## 716 Data availability statement 717 The APSIM wheat dataset analyzed in this study can be accessed via Bustos-Korts (2023). ## 718 Declaration of AI use 719 Parts of the **R** analysis script were written with the assistance of GitHub Copilot inline code completions 720 in RStudio 2024.12 during 2025. The completions were generated by GitHub’s default completion models 721 at the time: a GPT-3.5 Turbo-based model, GPT-4o Copilot and GPT-4.1 Copilot. The authors reviewed, 722 tested and revised all code suggested by Copilot and verified that the results are accurate and reproducible. 723 Generative AI was not used to design the model, interpret the results, or write the manuscript. AI 724 assistance (Claude Opus 5.5) was utilized during the final preparation of the submission strictly for 725 language enhancement, stylistic refinement, and copy editing. ## 726 Software 727 The computer code used in this paper is available through the dedicated GitHub repository at `[https:](https://github.com/gugushvili/dynamic_phenotypes)` 728 `[//github.com/gugushvili/dynamic_phenotypes](https://github.com/gugushvili/dynamic_phenotypes)` . ## 729 References 730 Adak, A., Kang, M., Anderson, S. L., Murray, S. C., Jarquin, D., Wong, R. K., and Katzfuß, M. Phenomic 731 data-driven biological prediction of maize through field-based high-throughput phenotyping integration 732 with genomic data. _Journal_ _of_ _Experimental_ _Botany_, 74(17):5307–5326, 2023. doi: 10.1093/jxb/erad216. 733 Araus, J. L., Kefauver, S. C., Zaman-Allah, M., Olsen, M. S., and Cairns, J. E. Translating high-throughput 734 phenotyping into genetic gain. _Trends_ _in_ _Plant_ _Science_, 23(5):451–466, 2018. doi: 10.1016/j.tplants. 735 2018.02.001. 736 Araus, J. L., Kefauver, S. C., Vergara-D´ıaz, O., Gracia-Romero, A., Rezzouk, F. Z., Segarra, J., Buchaillot, 737 M. L., Chang-Espino, M., Vatter, T., Sanchez-Bragado, R., et al. Crop phenotyping in a context of 738 global change: What to measure and how to do it. _Journal_ _of_ _Integrative_ _Plant_ _Biology_, 64(2):592–618, 739 2022. doi: 10.1111/jipb.13191. 740 Avagyan, V., Boer, M. P., Solin, J., van Dijk, A. D. J., Bustos-Korts, D., van Rossum, B.-J., Ramakers, J. 741 J. C., van Eeuwijk, F., and Kruijer, W. Penalized factorial regression as a flexible and computationally 742 attractive reaction norm model for prediction in the presence of _g × e_ . _Theoretical_ _and_ _Applied_ _Genetics_, 743 138(4):88, 2025. doi: 10.1007/s00122-025-04865-4. 744 Baba, T., Momen, M., Campbell, M. T., Walia, H., and Morota, G. Multi-trait random regression models 745 increase genomic prediction accuracy for a temporal physiological trait derived from high-throughput 746 phenotyping. _PLOS_ _ONE_, 15(2):e0228118, 2020. doi: 10.1371/journal.pone.0228118. 747 Baker, R. L., Leong, W. F., An, N., Brock, M. T., Rubin, M. J., Welch, S., and et al. Bayesian estimation 748 and use of high-throughput remote sensing indices for quantitative genetic analyses of leaf growth. 749 _Theor._ _Appl._ _Genet._, 131:283–298, 2018a. doi: 10.1007/s00122-017-3001-6. 750 Baker, R. L., Leong, W. F., Welch, S., and Weinig, C. Mapping and predicting non-linear Brassica rapa 751 growth phenotypes based on Bayesian and frequentist complex trait estimation. _G3:_ _Genes,_ _Genomes,_ 752 _Genetics_, 8:1247–1258, 2018b. doi: 10.1534/g3.117.300350. 17 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.05.764527; this version posted November 6, 2026. The copyright holder for this](https://doi.org/10.64898/2026.11.05.764527) preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 753 Baker, R. L. and Wang, D. R. Working with longitudinal data: quantifying developmental processes using 754 function-valued trait modeling. _American_ _Journal_ _of_ _Botany_, 108(6):905–908, 2021. doi: 10.1002/ajb2. 755 1677. 756 Bennett, S. Log-logistic regression models for survival data. _Journal_ _of_ _the_ _Royal_ _Statistical_ _Society._ 757 _Series_ _C_ _(Applied_ _Statistics)_, 32(2):165–171, 1983. ISSN 00359254, 14679876. URL `[http://www.jstor.](http://www.jstor.org/stable/2347295)` 758 `[org/stable/2347295](http://www.jstor.org/stable/2347295)` . 759 Betancourt, M. A conceptual introduction to Hamiltonian Monte Carlo. _arXiv_ _preprint_ _arXiv:1701.02434_, 760 2017. doi: 10.48550/arXiv.1701.02434. 761 Betancourt, M., Byrne, S., Livingstone, S., and Girolami, M. The geometric foundations of Hamiltonian 762 Monte Carlo. _Bernoulli_, 23(4A):2257–2298, 2017. doi: 10.3150/16-BEJ810. 763 Blei, D. M., Kucukelbir, A., and McAuliffe, J. D. Variational inference: a review for statisticians. _Journal_ 764 _of_ _the_ _American_ _Statistical_ _Association_, 112(518):859–877, 2017. doi: 10.1080/01621459.2017.1285773. 765 Bois, F. Y., Gelman, A., Jiang, J., Maszle, D. R., Zeise, L., and Alexeef, G. Population toxicokinetics of 766 tetrachloroethylene. _Archives_ _of_ _Toxicology_, 70(6):347–355, 1996. doi: 10.1007/s002040050284. 767 Bouidghaghen, J., Moreau, L., Beauchˆene, K., Chapuis, R., Mangel, N., Cabrera-Bosquet, L., Welcker, 768 C., Bogard, M., and Tardieu, F. Robotized indoor phenotyping allows genomic prediction of adaptive 769 traits in the field. _Nature_ _Communications_, 14(1):6603, 2023. doi: 10.1038/s41467-023-42298-z. 770 Brooks, S., Gelman, A., Jones, G. L., and Meng, X.-L., editors. _Handbook_ _of_ _Markov_ _chain_ _Monte_ 771 _Carlo_ . Chapman Hall/CRC Handb. Mod. Stat. Methods. Boca Raton, FL: CRC Press, 2011. ISBN 772 978-1-4200-7941-8; 978-1-4200-7942-5. doi: 10.1201/b10905. 773 B¨urkner, P.-C. brms: an R package for Bayesian multilevel models using Stan. _Journal_ _of_ _Statistical_ 774 _Software_, 80(1):1–28, 2017. doi: 10.18637/jss.v080.i01. 775 B¨urkner, P.-C. Advanced Bayesian multilevel modeling with the R package brms. _The_ _R_ _Journal_, 10(1): 776 395–411, 2018. doi: 10.32614/RJ-2018-017. 777 B¨urkner, P.-C. Bayesian item response modeling in R with brms and Stan. _Journal_ _of_ _Statistical_ _Software_, 778 100(5):1–54, 2021. doi: 10.18637/jss.v100.i05. 779 Bustos-Korts, D. APSIM-simulated wheat genotypes corresponding to the paper ‘Bustos-Korts, D. et al 780 2019. “From QTLs to Adaptation Landscapes: Using Genotype-To-Phenotype Models to Characterize 781 _G × E_ Over Time.” Frontiers in Plant Science 10:1–23. https://doi.org/10.3389/fpls.2019.01540.’. 12 782 2023. doi: 10.6084/m9.figshare.24843705.v1. 783 Bustos-Korts, D., Boer, M. P., Malosetti, M., Chapman, S., Chenu, K., Zheng, B., and van Eeuwijk, F. A. 784 Combining crop growth modeling and statistical genetic modeling to evaluate phenotyping strategies. 785 _Frontiers_ _in_ _Plant_ _Science_, 10:1491, 2019a. ISSN 1664-462X. doi: 10.3389/fpls.2019.01491. URL 786 `[https://www.frontiersin.org/articles/10.3389/fpls.2019.01491](https://www.frontiersin.org/articles/10.3389/fpls.2019.01491)` . 787 Bustos-Korts, D., Malosetti, M., Chenu, K., Chapman, S., Boer, M. P., Zheng, B., and van Eeuwijk, F. A. 788 From QTLs to adaptation landscapes: using genotype-to-phenotype models to characterize _G × E_ over 789 time. _Frontiers_ _in_ _Plant_ _Science_, 10:1540, 2019b. ISSN 1664-462X. doi: 10.3389/fpls.2019.01540. URL 790 `[https://www.frontiersin.org/articles/10.3389/fpls.2019.01540](https://www.frontiersin.org/articles/10.3389/fpls.2019.01540)` . 791 Calderhead, B. and Girolami, M. Estimating Bayes factors via thermodynamic integration and population 792 MCMC. _Comput._ _Stat._ _Data_ _Anal._, 53(12):4028–4045, 2009. ISSN 0167-9473. doi: 10.1016/j.csda.2009. 793 07.025. 794 Campbell, M., Walia, H., and Morota, G. Utilizing random regression models for genomic prediction of 795 a longitudinal trait derived from high-throughput phenotyping. _Plant_ _Direct_, 2(9):e00080, 2018. doi: 796 10.1002/pld3.80. 797 Casadebaig, P., Zheng, B., Chapman, S., Huth, N., Faivre, R., and Chenu, K. Assessment of the potential 798 impacts of wheat plant traits across environments by combining crop modeling and global sensitivity 799 analysis. _PLOS_ _ONE_, 11(1):1–27, 2016. doi: 10.1371/journal.pone.0146385. 18 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.05.764527; this version posted November 6, 2026. The copyright holder for this](https://doi.org/10.64898/2026.11.05.764527) preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 800 Chatzilena, A., van Leeuwen, E., Ratmann, O., Baguelin, M., and Demiris, N. Contemporary statistical 801 inference for infectious disease models using Stan. _Epidemics_, 29:100367, 2019. ISSN 1755-4365. 802 doi: 10.1016/j.epidem.2019.100367. URL `[https://www.sciencedirect.com/science/article/pii/](https://www.sciencedirect.com/science/article/pii/S1755436519300325)` 803 `[S1755436519300325](https://www.sciencedirect.com/science/article/pii/S1755436519300325)` . 804 Chen, Q., Zheng, B., Chen, T., and Chapman, S. C. Integrating a crop growth model and radiative 805 transfer model to improve estimation of crop traits based on deep learning. _Journal_ _of_ _Experimental_ 806 _Botany_, 73(19):6558–6574, 06 2022. ISSN 0022-0957. doi: 10.1093/jxb/erac291. 807 Chen, X., de S. N´oia-J´unior, R., Raihauti, T., Manceau, L., Dueri, S., Ramakers, J. J., Malik, W. A., 808 Bustos-Korts, D., Flamm, C., Grizeau, C., Herrera, J. M., Hiltbrunner, J., Hor´akov´a’, V., Levy 809 H¨aner, L., Masson, F., P´ecs, M., Povoln´y, M., Starnberger, P., Visse-Mansiaux, M., Collonnier, C., 810 Laurens, F., Piepho, H.-P., van Eeuwijk, F. A., Welcker, C., Parent, B., and Martre, P. Envirotyping 811 framework for wheat in Europe. _Field_ _Crops_ _Research_, 346:110616, 2026. ISSN 0378-4290. doi: 812 10.1016/j.fcr.2026.110616. 813 Chenu, K., Deihimfard, R., and Chapman, S. C. Large-scale characterization of drought pattern: a 814 continent-wide modelling approach applied to the Australian wheatbelt — spatial and temporal trends. 815 _New_ _Phytologist_, 198(3):801–820, 2013. doi: 10.1111/nph.12192. URL `[https://nph.onlinelibrary.](https://nph.onlinelibrary.wiley.com/doi/abs/10.1111/nph.12192)` 816 `[wiley.com/doi/abs/10.1111/nph.12192](https://nph.onlinelibrary.wiley.com/doi/abs/10.1111/nph.12192)` . 817 Chopin, N. and Papaspiliopoulos, O. _An_ _introduction_ _to_ _sequential_ _Monte_ _Carlo_ . Springer Ser. Stat. 818 Cham: Springer, 2020. ISBN 978-3-030-47844-5; 978-3-030-47847-6; 978-3-030-47845-2. doi: 10.1007/ 819 978-3-030-47845-2. 820 Comets, E., Lavenu, A., and Lavielle, M. Parameter estimation in nonlinear mixed effect models using 821 saemix, an R implementation of the SAEM algorithm. _Journal_ _of_ _Statistical_ _Software_, 80(3):1–41, 822 2017. doi: 10.18637/jss.v080.i03. URL `[https://www.jstatsoft.org/index.php/jss/article/view/](https://www.jstatsoft.org/index.php/jss/article/view/v080i03)` 823 `[v080i03](https://www.jstatsoft.org/index.php/jss/article/view/v080i03)` . 824 Cooper, M., Technow, F., Messina, C., Gho, C., and Totir, L. R. Use of crop growth models with 825 whole-genome prediction: application to a maize multienvironment trial. _Crop_ _Science_, 56(5):2141–2156, 826 2016. doi: 10.2135/cropsci2015.08.0512. 827 De Jong, P. and Shephard, N. The simulation smoother for time series models. _Biometrika_, 82(2):339–350, 828 1995. ISSN 0006-3444. doi: 10.1093/biomet/82.2.339. 829 Donnet, S. and Samson, A. Estimation of parameters in incomplete data models defined by dynam830 ical systems. _Journal_ _of_ _Statistical_ _Planning_ _and_ _Inference_, 137(9):2815–2831, 2007. ISSN 0378831 3758. doi: 10.1016/j.jspi.2006.10.013. URL `[https://www.sciencedirect.com/science/article/](https://www.sciencedirect.com/science/article/pii/S0378375807000353)` 832 `[pii/S0378375807000353](https://www.sciencedirect.com/science/article/pii/S0378375807000353)` . 833 Driver, C. C. and Voelkle, M. C. Hierarchical Bayesian continuous time dynamic modeling. _Psychological_ 834 _Methods_, 23(4):774–799, 2018. doi: 10.1037/met0000168. 835 Faraway, J. J. _Extending_ _the_ _linear_ _model_ _with_ _R._ _Generalized_ _linear,_ _mixed_ _effects_ _and_ _nonparametric_ 836 _regression_ _models._ Chapman Hall/CRC Texts Stat. Sci. Ser. Boca Raton, FL: CRC Press, 2nd 837 edition edition, 2016. ISBN 978-1-4987-2096-0; 978-1-4987-2098-4. doi: 10.1201/9781315382722. URL 838 `www.taylorfrancis.com/books/9781498720984` . 839 Fernandez, J. A., Hu, P., Arief, V. N., Gho, C., Zheng, B., Ramakers, J. J., Boer, M. P., van Eeuwijk, F. A., 840 and Chapman, S. C. Envirotypes for wheat _G × E_ in Australia: integrating multiple environmental 841 indices to characterize variety trials. _G3_ _Genes—Genomes—Genetics_, page jkag222, 08 2026. ISSN 842 2160-1836. doi: 10.1093/g3journal/jkag222. 843 Focke, W. W., van der Westhuizen, I., Musee, N., and Loots, M. T. Kinetic interpretation of log-logistic 844 dose-time response curves. _Scientific_ _Reports_, 7(1):2234, 2017. doi: 10.1038/s41598-017-02474-w. 845 Gabry, J., Mahr, T., B¨urkner, P.-C., Modr´ak, M., Barrett, M., Weber, F., Coronado Sroka, E., Sailynoja, 846 T., and Vehtari, A. Visual MCMC diagnostics. _URL_ _https://CRAN.R-Project.org/package=bayesplot._ 847 _Vignette_ _included_ _in_ _R_ _package_ _bayesplot,_ _version_ _1.13.0_, 2025. 848 Gelfand, A. E. and Smith, A. F. M. Sampling-based approaches to calculating marginal densities. _J._ _Am._ 849 _Stat._ _Assoc._, 85(410):398–409, 1990. ISSN 0162-1459. doi: 10.2307/2289776. 19 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.05.764527; this version posted November 6, 2026. The copyright holder for this](https://doi.org/10.64898/2026.11.05.764527) preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 850 Gelman, A. and Hill, J. _Data_ _analysis_ _using_ _regression_ _and_ _multilevel/hierarchical_ _models_ . Analytical 851 Methods for Social Research. Cambridge University Press, 2006. doi: 10.1017/CBO9780511790942. 852 Gelman, A. and Rubin, D. B. Inference from iterative simulation using multiple sequences. _Stat._ _Sci._, 7 853 (4):457–472, 1992. ISSN 0883-4237. doi: 10.1214/ss/1177011136. 854 Gelman, A., Bois, F., and Jiang, J. Physiological pharmacokinetic analysis using population modeling 855 and informative prior distributions. _J._ _Am._ _Stat._ _Assoc._, 91(436):1400–1412, 1996. ISSN 0162-1459. doi: 856 10.2307/2291566. 857 Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., and Rubin, D. B. _Bayesian_ _data_ 858 _analysis_ . Chapman Hall/CRC Texts Stat. Sci. Ser. Boca Raton, FL: CRC Press, 3rd revised and 859 updated edition, 2014. ISBN 978-1-4398-4095-5. doi: 10.1201/b16018. 860 Gelman, A., Simpson, D., and Betancourt, M. The prior can often only be understood in the context 861 of the likelihood. _Entropy_, 19(10), 2017. ISSN 1099-4300. doi: 10.3390/e19100555. URL `[https:](https://www.mdpi.com/1099-4300/19/10/555)` 862 `[//www.mdpi.com/1099-4300/19/10/555](https://www.mdpi.com/1099-4300/19/10/555)` . 863 Gelman, A., Hill, J., and Vehtari, A. _Regression_ _and_ _other_ _stories_ . Anal. Methods Soc. Res. Cambridge: 864 Cambridge University Press, 2021. ISBN 978-1-107-67651-0; 978-1-107-02398-7; 978-1-139-16187-9. doi: 865 10.1017/9781139161879. 866 Geman, S. and Geman, D. Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of 867 images. _IEEE_ _Trans._ _Pattern_ _Anal._ _Mach._ _Intell._, 6:721–741, 1984. ISSN 0162-8828. doi: 10.1109/ 868 TPAMI.1984.4767596. 869 Girolami, M. Bayesian inference for differential equations. _Theor._ _Comput._ _Sci._, 408(1):4–16, 2008. ISSN 870 0304-3975. doi: 10.1016/j.tcs.2008.07.005. 871 G´omez-Rubio, V. _Bayesian_ _inference_ _with_ _INLA_ . Boca Raton, FL: CRC Press, 2020. ISBN 978-1-138872 03987-2; 978-1-032-17453-2; 978-1-315-17558-4. doi: 10.1201/9781315175584. 873 Gomulkiewicz, R., Kingsolver, J. G., Carter, P. A., and Heckman, N. Variation and evolution of 874 function-valued traits. _Annual_ _Review_ _of_ _Ecology,_ _Evolution,_ _and_ _Systematics_, 49(1):139–164, 2018. doi: 875 10.1146/annurev-ecolsys-110316-022830. 876 Grinsztajn, L., Semenova, E., Margossian, C. C., and Riou, J. Bayesian workflow for disease transmission 877 modeling in Stan. _Statistics_ _in_ _Medicine_, 40(27):6209–6234, 2021. doi: 10.1002/sim.9164. 878 Gu´edon, Y., Caraglio, Y., Granier, C., Lauri, P.-E., [´] and Muller, B. Identifying developmental patterns 879 in structured plant phenotyping data. _Plant_ _Systems_ _Biology:_ _Methods_ _and_ _Protocols_, pages 199–225, 880 2022. doi: 10.1007/978-1-0716-1816-5 ~~1~~ 0. 881 Hoffman, M. D. and Gelman, A. The No-U-Turn Sampler: adaptively setting path lengths in Hamiltonian 882 Monte Carlo. _Journal_ _of_ _Machine_ _Learning_ _Research_, 15(47):1593–1623, 2014. URL `[http://jmlr.org/](http://jmlr.org/papers/v15/hoffman14a.html)` 883 `[papers/v15/hoffman14a.html](http://jmlr.org/papers/v15/hoffman14a.html)` . 884 Holzworth, D., Huth, N. I., Fainges, J., Brown, H., Zurcher, E., Cichota, R., Verrall, S., Herrmann, 885 N. I., Zheng, B., and Snow, V. APSIM next generation: overcoming challenges in modernising a 886 farming systems model. _Environmental_ _Modelling_ _&_ _Software_, 103:43–51, 2018. ISSN 1364-8152. 887 doi: 10.1016/j.envsoft.2018.02.002. URL `[https://www.sciencedirect.com/science/article/pii/](https://www.sciencedirect.com/science/article/pii/S1364815217311921)` 888 `[S1364815217311921](https://www.sciencedirect.com/science/article/pii/S1364815217311921)` . 889 Jarqu´ın, D., Crossa, J., Lacaze, X., Du Cheyron, P., Daucourt, J., Lorgeou, J., Piraux, F., Guerreiro, 890 L., P´erez, P., Calus, M., et al. A reaction norm model for genomic selection using high-dimensional 891 genomic and environmental data. _Theoretical_ _and_ _Applied_ _Genetics_, 127:595–607, 2014. doi: 10.1007/ 892 s00122-013-2243-1. 893 Jighly, A., Weeks, A., Christy, B., O’Leary, G. J., Kant, S., Aggarwal, R., Hessel, D., Forrest, K. L., 894 Technow, F., Tibbits, J. F. G., Totir, R., Spangenberg, G. C., Hayden, M. J., Munkvold, J., and 895 Daetwyler, H. D. Integrating biophysical crop growth models and whole genome prediction for their 896 mutual benefit: a case study in wheat phenology. _Journal_ _of_ _Experimental_ _Botany_, 74(15):4415–4426, 897 05 2023. ISSN 0022-0957. doi: 10.1093/jxb/erad162. 20 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.05.764527; this version posted November 6, 2026. The copyright holder for this](https://doi.org/10.64898/2026.11.05.764527) preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 898 Kaheman, K., Kutz, J. N., and Brunton, S. L. Sindy-pi: a robust algorithm for parallel implicit sparse 899 identification of nonlinear dynamics. _Proceedings_ _of_ _the_ _Royal_ _Society_ _A_, 476(2242):20200279, 2020. doi: 900 10.1098/rspa.2020.0279. 901 Kay, M. tidybayes: tidy data and geoms for Bayesian models. _URL_ _http://mjskay.github.io/tidybayes._ 902 _R_ _package_ _version_ _3.0.2_, 2022. doi: 10.5281/zenodo.1308151. URL `[http://mjskay.github.io/](http://mjskay.github.io/tidybayes/)` 903 `[tidybayes/](http://mjskay.github.io/tidybayes/)` . 904 Keating, B. A., Carberry, P. S., Hammer, G. L., Probert, M. E., Robertson, M. J., Holzworth, D., Huth, 905 N. I., Hargreaves, J. N. G., Meinke, H., Hochman, Z., McLean, G., Verburg, K., Snow, V., Dimes, J. P., 906 Silburn, M., Wang, E., Brown, S., Bristow, K. L., Asseng, S., Chapman, S., McCown, R. L., Freebairn, 907 D. M., and Smith, C. J. An overview of APSIM, a model designed for farming systems simulation. 908 _European_ _Journal_ _of_ _Agronomy_, 18(3):267–288, 2003. ISSN 1161-0301. doi: 10.1016/S1161-0301(02) 909 00108-9. URL `[https://www.sciencedirect.com/science/article/pii/S1161030102001089](https://www.sciencedirect.com/science/article/pii/S1161030102001089)` . 910 Kruschke, J. _Doing_ _Bayesian_ _data_ _analysis._ _A_ _tutorial_ _introduction_ _with_ _R,_ _JAGS,_ _and_ _Stan_ . Amsterdam: 911 Elsevier/Academic Press, 2nd edition, 2014. ISBN 978-0-12-405888-0; 978-0-12-405916-0. URL `www.` 912 `sciencedirect.com/science/book/9780124058880` . 913 Lee, S. Y. Bayesian nonlinear models for repeated measurement data: An overview, implementation, 914 and applications. _Mathematics_, 10(6), 2022. ISSN 2227-7390. doi: 10.3390/math10060898. URL 915 `[https://www.mdpi.com/2227-7390/10/6/898](https://www.mdpi.com/2227-7390/10/6/898)` . 916 Li, Z. and Sillanp¨a¨a, M. J. Overview of LASSO-related penalized regression methods for quantitative 917 trait mapping and genomic selection. _Theoretical_ _and_ _Applied_ _Genetics_, 125:419–435, 2012. doi: 918 10.1007/s00122-012-1892-9. 919 Lindsten, F. and Sch¨on, T. B. Backward simulation methods for Monte Carlo statistical inference. _Found._ 920 _Trends_ _Mach._ _Learn._, 6(1):1–143, 2013. ISSN 1935-8237. doi: 10.1561/2200000045. 921 MacQueen, J. Some methods for classification and analysis of multivariate observations. In Le Cam, L. M. 922 and Neyman, J., editors, _Proc._ _5th_ _Berkeley_ _Symp._ _Math._ _Stat._ _Probab.,_ _Univ._ _Calif_ _(Berkeley,_ _Calif.,_ 923 _1965/66)_, volume 1, pages 281–297. University of California Press, 1967. 924 Messina, C. D., Technow, F., Tang, T., Totir, R., Gho, C., and Cooper, M. Leveraging biological 925 insight and environmental variation to improve phenotypic prediction: Integrating crop growth models 926 (cgm) with whole genome prediction (wgp). _European_ _Journal_ _of_ _Agronomy_, 100:151–162, 2018. doi: 927 10.1016/j.eja.2018.01.007. 928 Momen, M., Campbell, M. T., Walia, H., and Morota, G. Predicting longitudinal traits derived from 929 high-throughput phenomics in contrasting environments using genomic Legendre polynomials and 930 B-splines. _G3:_ _Genes,_ _Genomes,_ _Genetics_, 9(10):3369–3380, 2019. doi: 10.1534/g3.119.400346. 931 Ogle, K. and Barber, J. J. Ensuring identifiability in hierarchical mixed effects Bayesian models. _Ecological_ 932 _Applications_, 30(7):e02159, 2020. doi: 10.1002/eap.2159. 933 Onogi, A. Connecting mathematical models to genomes: joint estimation of model parameters and 934 genome-wide marker effects on these parameters. _Bioinformatics_, 36(10):3169–3176, 02 2020. ISSN 935 1367-4803. doi: 10.1093/bioinformatics/btaa129. 936 Onogi, A., Ideta, O., Inoshita, Y., Ebana, K., Yoshioka, T., Yamasaki, M., and Iwata, H. Exploring the 937 areas of applicability of whole-genome prediction methods for asian rice (oryza sativa l.). _Theoretical_ 938 _and_ _Applied_ _Genetics_, 128(1):41–53, 2015. doi: 10.1007/s00122-014-2411-y. 939 Onogi, A., Ogino, A., Sato, A., Kurogi, K., Yasumori, T., and Togashi, K. Development of a structural 940 growth curve model that considers the causal effect of initial phenotypes. _Genetics_ _Selection_ _Evolution_, 941 51(1):19, 2019. doi: 10.1186/s12711-019-0461-y. 942 Panday, U. S., Pratihast, A. K., Aryal, J., and Kayastha, R. B. A review on drone-based data solutions 943 for cereal crops. _Drones_, 4(3), 2020a. ISSN 2504-446X. doi: 10.3390/drones4030041. URL `[https:](https://www.mdpi.com/2504-446X/4/3/41)` 944 `[//www.mdpi.com/2504-446X/4/3/41](https://www.mdpi.com/2504-446X/4/3/41)` . 21 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.05.764527; this version posted November 6, 2026. The copyright holder for this](https://doi.org/10.64898/2026.11.05.764527) preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 945 Panday, U. S., Shrestha, N., Maharjan, S., Pratihast, A. K., Shahnawaz, Shrestha, K. L., and Aryal, 946 J. Correlating the plant height of wheat with above-ground biomass and crop yield using drone 947 imagery and crop surface model, a case study from Nepal. _Drones_, 4(3), 2020b. ISSN 2504-446X. doi: 948 10.3390/drones4030028. URL `[https://www.mdpi.com/2504-446X/4/3/28](https://www.mdpi.com/2504-446X/4/3/28)` . 949 P´erez-Valencia, D. M., Rodr´ıguez-Alvarez, M. X., Boer, M. P., Kronenberg, L., Hund, A., Cabrera-Bosquet, [´] 950 L., Millet, E. J., and Eeuwijk, F. A. v. A two-stage approach for the spatio-temporal analysis of high951 throughput phenotyping data. _Scientific_ _Reports_, 12(1):3177, 2022. doi: 10.1038/s41598-022-06935-9. 952 P´erez-Valencia, D. M., Rodr´ıguez-Alvarez, [´] M. X., Boer, M. P., and van Eeuwijk, F. A. A one-stage 953 approach for the spatio-temporal analysis of high-throughput phenotyping data. _Journal_ _of_ _Agricultural,_ 954 _Biological_ _and_ _Environmental_ _Statistics_, pages 1–23, 2024. doi: 10.1007/s13253-024-00642-w. 955 Pinheiro, J. C. and Bates, D. M. _Mixed-effects_ _models_ _in_ _S_ _and_ _S-Plus_ . Stat. Comput. (Cham). New York, 956 NY: Springer, 2000. ISBN 0-387-98957-9. doi: 10.1007/b98882. 957 Plan, E. L., Maloney, A., Mentr´e, F., Karlsson, M. O., and Bertrand, J. Performance comparison of 958 various maximum likelihood nonlinear mixed-effects estimation methods for dose-response models. _The_ 959 _AAPS_ _Journal_, 14(3):420–432, 2012. doi: 10.1208/s12248-012-9349-2. 960 Poudel, P., Naidenov, B., Chen, C., Alderman, P. D., and Welch, S. M. Integrating genomic prediction 961 and genotype specific parameter estimation in ecophysiological models: overview and perspectives. _in_ 962 _silico_ _Plants_, 5(1):diad007, 2023. doi: 10.1093/insilicoplants/diad007. 963 R Core Team. _R:_ _a_ _language_ _and_ _environment_ _for_ _statistical_ _computing_ . R Foundation for Statistical 964 Computing, Vienna, Austria, 2022. URL `[https://www.R-project.org/](https://www.R-project.org/)` . 965 Ramsay, J. and Hooker, G. _Dynamic_ _data_ _analysis._ _Modeling_ _data_ _with_ _differential_ _equations_ . Springer 966 Ser. Stat. New York, NY: Springer, 2017. ISBN 978-1-4939-7188-6; 978-1-4939-7190-9. doi: 10.1007/ 967 978-1-4939-7190-9. 968 Resende, R. T., Piepho, H.-P., Rosa, G. J., Silva-Junior, O. B., e Silva, F. F., de Resende, M. D. V., and 969 Grattapaglia, D. Enviromics in breeding: applications and perspectives on envirotypic-assisted selection. 970 _Theoretical_ _and_ _Applied_ _Genetics_, 134:95–112, 2021. doi: 10.1007/s00122-020-03684-z. 971 Ritz, C., Jensen, S. M., Gerhard, D., and Streibig, J. C. _Dose-response_ _analysis_ _using_ _R_ . Chapman 972 Hall/CRC The R Ser. Boca Raton, FL: CRC Press, 2020. ISBN 978-1-138-03431-0; 978-1-315-27009-8. 973 doi: 10.1201/b21966. 974 Robert, P., Goudemand, E., Auzanneau, J., Oury, F.-X., Rolland, B., Heumez, E., Bouchet, S., Caillebotte, 975 A., Mary-Huard, T., Le Gouis, J., and Rincent, R. Phenomic selection in wheat breeding: prediction of 976 the genotype-by-environment interaction in multi-environment breeding trials. _Theoretical_ _and_ _Applied_ 977 _Genetics_, 135(10):3337–3356, 2022. doi: 10.1007/s00122-022-04170-4. 978 Roberts, G. O. and Rosenthal, J. S. General state space Markov chains and MCMC algorithms. _Probab._ 979 _Surv._, 1:20–71, 2004. ISSN 1549-5787. doi: 10.1214/154957804100000024. URL `[https://eudml.org/](https://eudml.org/doc/224487)` 980 `[doc/224487](https://eudml.org/doc/224487)` . 981 Ross, G. J. S. _Nonlinear_ _estimation_ . Springer Ser. Stat. New York etc.: Springer-Verlag, 1990. ISBN 982 3-540-97278-1. doi: 10.1007/978-1-4612-3412-8. 983 Roth, L., Rodr´ıguez-Alvarez, [´] M. X., van Eeuwijk, F., Piepho, H.-P., and Hund, A. Phenomics data 984 processing: A plot-level model for repeated measurements to extract the timing of key stages and 985 quantities at defined time points. _Field_ _Crops_ _Research_, 274:108314, 2021. doi: 10.1016/j.fcr.2021.108314. 986 Roth, L., Piepho, H.-P., and Hund, A. Phenomics data processing: extracting dose–response curve param987 eters from high-resolution temperature courses and repeated field-based wheat height measurements. _in_ 988 _silico_ _Plants_, 4(1):diac007, 2022. doi: 10.1093/insilicoplants/diac007. 989 Roth, L., Fossati, D., Kr¨ahenb¨uhl, P., Walter, A., and Hund, A. Image-based phenomic prediction can 990 provide valuable decision support in wheat breeding. _Theoretical_ _and_ _Applied_ _Genetics_, 136(7):162, 991 2023. doi: 10.1007/s00122-023-04395-x. 992 Searle, S. R., Casella, G., and McCulloch, C. E. _Variance_ _components_ . Wiley Ser. Probab. Stat. Hoboken, 993 NJ: John Wiley & Sons, 2006. ISBN 0-470-00959-4. doi: 10.1002/9780470316856. 22 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.05.764527; this version posted November 6, 2026. The copyright holder for this](https://doi.org/10.64898/2026.11.05.764527) preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 994 Seber, G. A. F. and Wild, C. J. _Nonlinear_ _regression_ . Wiley Ser. Probab. Math. Stat. Chichester etc.: 995 John Wiley & Sons Ltd., 1989. ISBN 0-471-61760-1. doi: 10.1002/0471725315. 996 Shao, Y., van Eeuwijk, F., Peeters, C. F., Zumsteg, O., Athanasiadis, I. N., and van Voorn, G. Physics997 Informed Neural Network methods for predicting plant height development. _Computers_ _and_ _Electronics_ 998 _in_ _Agriculture_, 251:111988, 2026. ISSN 0168-1699. doi: 10.1016/j.compag.2026.111988. 999 Sisson, S. A., Fan, Y., and Beaumont, M., editors. _Handbook_ _of_ _approximate_ _Bayesian_ _computation_ . 1000 Chapman Hall/CRC Handb. Mod. Stat. Methods. Boca Raton, FL: CRC Press, 2019. ISBN 978-11001 4398-8150-7; 978-1-315-11719-5. doi: 10.1201/9781315117195. 1002 Smith, D. T., Potgieter, A. B., and Chapman, S. C. Scaling up high-throughput phenotyping for 1003 abiotic stress selection in the field. _Theoretical_ _and_ _Applied_ _Genetics_, 134(6):1845–1866, 2021. doi: 1004 10.1007/s00122-021-03864-5. 1005 Stan Development Team. _Stan_ _modeling_ _language_ _user’s_ _guide_ _and_ _reference_ _manual._ _Version_ _2.33._ 1006 _https://mc-stan.org_ . NumFOCUS, 2023. URL `[https://mc-stan.org](https://mc-stan.org)` . ˇ 1007 Strumbelj, E., Bouchard-Cˆot´e, A., Corander, J., Gelman, A., Rue, H., Murray, L., Pesonen, H., Plummer, 1008 M., and Vehtari, A. Past, present, and future of software for Bayesian inference. _Statistical_ _Science_, 39 1009 (1):46–61, 2024. doi: 10.1214/23-STS907. 1010 Sun, J., Rutkoski, J. E., Poland, J. A., Crossa, J., Jannink, J.-L., and Sorrells, M. E. Multitrait, 1011 random regression, or simple repeatability model in high-throughput phenotyping data improve genomic 1012 prediction for wheat grain yield. _The_ _Plant_ _Genome_, 10(2):plantgenome2016.11.0111, 2017. doi: 1013 10.3835/plantgenome2016.11.0111. 1014 Sun, J., Poland, J. A., Mondal, S., Crossa, J., Juliana, P., Singh, R. P., Rutkoski, J. E., Jannink, J.-L., 1015 Crespo-Herrera, L., Velu, G., et al. High-throughput phenotyping platforms enhance genomic selection 1016 for wheat grain yield across populations and cycles in early stage. _Theoretical_ _and_ _Applied_ _Genetics_, 1017 132:1705–1720, 2019. doi: 10.1007/s00122-019-03309-0. 1018 Tardieu, F., Cabrera-Bosquet, L., Pridmore, T., and Bennett, M. Plant phenomics, from sensors to 1019 knowledge. _Current_ _Biology_, 27(15):R770–R783, 2017. doi: 10.1016/j.cub.2017.05.055. 1020 Technow, F., Messina, C. D., Totir, L. R., and Cooper, M. Integrating crop growth models with whole 1021 genome prediction through Approximate Bayesian Computation. _PLOS_ _ONE_, 10(6):1–20, 06 2015. doi: 1022 10.1371/journal.pone.0130855. 1023 Transtrum, M. K., Machta, B. B., and Sethna, J. P. Geometry of nonlinear least squares with applications to 1024 sloppy models and optimization. _Phys._ _Rev._ _E_, 83:036701, Mar 2011. doi: 10.1103/PhysRevE.83.036701. 1025 URL `[https://link.aps.org/doi/10.1103/PhysRevE.83.036701](https://link.aps.org/doi/10.1103/PhysRevE.83.036701)` . 1026 Van der Laan, M., Pollard, K., and Bryan, J. A new partitioning around medoids algorithm. _Journal_ _of_ 1027 _Statistical_ _Computation_ _and_ _Simulation_, 73(8):575–584, 2003. doi: 10.1080/0094965031000136012. 1028 Van Eeuwijk, F. A., Bustos-Korts, D. V., and Malosetti, M. What should students in plant breeding know 1029 about the statistical aspects of genotype _×_ environment interactions? _Crop_ _Science_, 56(5):2119–2140, 1030 2016. doi: 10.2135/cropsci2015.06.0375. 1031 Van Voorn, G. A. K., Boer, M. P., Truong, S. H., Friedenberg, N. A., Gugushvili, S., McCormick, R., 1032 Bustos Korts, D., Messina, C. D., and Van Eeuwijk, F. A. A conceptual framework for the dynamic 1033 modeling of time-resolved phenotypes for sets of genotype-environment-management combinations (I): 1034 a model library. _Front._ _Plant_ _Sci._, 14:1172359, 2023. ISSN 1664-462X. doi: 10.3389/fpls.2023.1172359. 1035 Varona, L., Moreno, C., Cort´es, L. A. G., Yag¨ue, G., and Altarriba, J. Two-step versus joint analysis 1036 of Von Bertalanffy function. _Journal_ _of_ _Animal_ _Breeding_ _and_ _Genetics_, 116(5):331–338, 1999. doi: 1037 10.1046/j.1439-0388.1999.00220.x. 1038 Vehtari, A., Gelman, A., Simpson, D., Carpenter, B., and B¨urkner, P.-C. Rank-normalization, folding, 1039 and localization: an improved _R_ [ˆ︁] for assessing convergence of MCMC (with Discussion). _Bayesian_ _Anal._, 1040 16(2):667–718, 2021. ISSN 1931-6690. doi: 10.1214/20-BA1221. 1041 Verburg, C., Heinlein, A., van Eeuwijk, F., and Budko, N. Physiology-Informed Sparse Identification of 1042 Dynamic Crop Models. _arXiv_ _e-prints_, art. arXiv:2610.18920, oct 2026. doi: 10.48550/arXiv.2610.18920. 23 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.05.764527; this version posted November 6, 2026. The copyright holder for this](https://doi.org/10.64898/2026.11.05.764527) preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 1043 Voit, E. O. and Almeida, J. Decoupling dynamical systems for pathway identification from metabolic 1044 profiles. _Bioinformatics_, 20(11):1670–1681, 02 2004. ISSN 1367-4803. doi: 10.1093/bioinformatics/ 1045 bth140. 1046 Wakefield, J. The Bayesian analysis of population pharmacokinetic models. _J._ _Am._ _Stat._ _Assoc._, 91(433): 1047 62–75, 1996. ISSN 0162-1459. doi: 10.2307/2291383. 1048 Wang, Z., Pang, X., Wu, W., Wang, J., Wang, Z., and Wu, R. Modeling phenotypic plasticity in growth 1049 trajectories: a statistical framework. _Evolution_, 68(1):81–91, 2014. doi: 10.1111/evo.12263. 1050 Wickham, H. _ggplot2:_ _elegant_ _graphics_ _for_ _data_ _analysis_ . Springer-Verlag New York, 2016. ISBN 1051 978-3-319-24277-4. URL `[https://ggplot2.tidyverse.org](https://ggplot2.tidyverse.org)` . 1052 Xu, Y. Envirotyping for deciphering environmental impacts on crop plants. _Theoretical_ _and_ _Applied_ 1053 _Genetics_, 129:653–673, 2016. doi: 10.1007/s00122-016-2691-5. 1054 Zheng, B., Biddulph, B., Li, D., Kuchel, H., and Chapman, S. Quantification of the effects of VRN1 and 1055 Ppd-D1 to predict spring wheat (Triticum aestivum) heading time across diverse environments. _Journal_ 1056 _of_ _Experimental_ _Botany_, 64(12):3747–3761, 07 2013. ISSN 0022-0957. doi: 10.1093/jxb/ert209. 24 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.05.764527; this version posted November 6, 2026. The copyright holder for this](https://doi.org/10.64898/2026.11.05.764527) preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. ## 1057 5 Figures Figure 1: The Kruschke style diagram for the hierarchical Bayesian model from Section 2.5. See the main text for the notation. The inner and outer layers actually consist of three identical copies of distributions, although only one is shown. For each coordinate _θ_ _ij_ _[p]_ [of] [the] [parameter] [vector] _[θ]_ _[ij]_ [,] [with] [the] [superscript] varying over _p_ = _M_, _mid_ and _r_, a separate distribution is specified. The indices halfway the arrows indicate the way in which the variables or parameters are replicated over those indices. 25 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.05.764527; this version posted November 6, 2026. The copyright holder for this](https://doi.org/10.64898/2026.11.05.764527) preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. Figure 2: Observed daily APSIM wheat biomass data for medoid year at each location showing all 48 genotypes. Lines are colored according to genotypes. Colors are consistent across different panels. The red tick marks along the horizontal axis indicate time points at which biomass observations were used for Bayesian hierarchical modeling. 26 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.05.764527; this version posted November 6, 2026. The copyright holder for this](https://doi.org/10.64898/2026.11.05.764527) preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. Figure 3: Fitted curves (posterior medians) for medoid year at each location showing all genotypes. Points represent observations. Lines and points are colored according to genotypes. Colors are consistent across panels. The time axis has been obtained by dropping the first 14 days after sowing and rescaling by division by 10. 27 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.05.764527; this version posted November 6, 2026. The copyright holder for this](https://doi.org/10.64898/2026.11.05.764527) preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. Figure 4: Fitted curves (posterior medians) for medoid genotype in medoid year for each of 4 locations. Observations at locations are distinguished by symbol: points, squares, triangles and crosses. Colored bands represent 95% Bayesian (pointwise) prediction bands. Bands are colored according to locations. The time axis has been obtained by dropping the first 14 days after sowing and rescaling by division by 10. 28 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.05.764527; this version posted November 6, 2026. The copyright holder for this](https://doi.org/10.64898/2026.11.05.764527) preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. Figure 5: Fitted curves (posterior medians) for the medoid genotype across all years at each of 4 locations. Points represent observations. Lines and points are colored according to years. Colors are consistent across panels. The time axis has been obtained by dropping the first 14 days after sowing and rescaling by division by 10. 29 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.05.764527; this version posted November 6, 2026. The copyright holder for this](https://doi.org/10.64898/2026.11.05.764527) preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. _α_ _[p]_ _[,]_ _[σ]_ _[p]_ Figure 6: Posterior summaries for parameters _σ_ _α_ _[p]_ _[,]_ _[σ]_ _β_ _[p]_ [,] [and] _[σ]_ ( _[p]_ _α·β_ ) [for] [each] [of] [the] [superscripts] _[p]_ [ =] _[ M]_ [,] _[mid]_ and _r_ . Each panel gathers results for a different parameter of the log-logistic function, identified by labels `Asymptote`, `Midpoint` and `Slope` . Points in the middle of the segments are posterior medians. The thick and thin segments represent 67% and 95% credible intervals, respectively. No posterior densities are plotted, as they are highly peaked for _σ_ _α_ _[p]_ [and] _[σ]_ _[p]_ _·_ [parameters,] [and] [as] [such] [visually] [uninformative.] Figure 6: Posterior summaries for parameters _σ_ _α_ _[p]_ _[p]_ [and] _[σ]_ _[p]_ _β_ [,] ( _α_ _[p]_ [and] _[σ]_ _[p]_ ( _[p]_ _α·β_ ) [parameters,] [and] [as] [such] [visually] [uninformative.] 30 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.05.764527; this version posted November 6, 2026. The copyright holder for this](https://doi.org/10.64898/2026.11.05.764527) preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. Figure 7: Posterior summaries for ratios _σ_ _[p]_ ( _[p]_ _α·β_ ) _[/σ]_ _α_ _[p]_ Figure 7: Posterior summaries for ratios _σ_ ( _[p]_ _α·β_ ) _[/σ]_ _α_ _[p]_ [for] [each] [of] [the] [superscripts] _[p]_ [ =] _[ M]_ [,] _[mid]_ [and] _[r]_ [.] [Labels] `Asymptote`, `Midpoint` and `Slope` refer to the parameters of the log-logistic function. Points in the middle of the segments are posterior medians. The thick and thin segments represent 67% and 95% credible intervals, respectively. The densities are the posterior densities. 31 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.05.764527; this version posted November 6, 2026. The copyright holder for this](https://doi.org/10.64898/2026.11.05.764527) preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. Figure 8: Percentages of genotypes in individual trials with interactions ( _α · β_ ) _[M]_ _ij_ [(for] [the] [asymptote] parameter _M_ _[max]_ ) whose 95% credible interval excludes zero. The higher this percentage, the more specific are the environmental challenges imposed by the conditions. 32 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.05.764527; this version posted November 6, 2026. The copyright holder for this](https://doi.org/10.64898/2026.11.05.764527) preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. Figure 9: Percentage of interactions for individual genotypes ( _α · β_ ) _[M]_ _ij_ [across] [all] [trials,] [whose] [95%] [credible] intervals exclude zero. This percentage can be read as a measure of stability. The lower the percentage, the more stable is a genotype. 33 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.05.764527; this version posted November 6, 2026. The copyright holder for this](https://doi.org/10.64898/2026.11.05.764527) preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. Figure 10: Correlations between APSIM parameters, on the one hand, and the genotype main effects _α_ _i_ _[p]_ [,] [on] [the] [other,] [for] [superscripts] _[p]_ [=] _[M]_ [,] _[mid]_ [and] _[r]_ [.] [APSIM] [parameter] [labels] [are] [explained] [in] [Table] 4. Each panel gathers results for a different parameter of the log-logistic function, identified by labels `Asymptote`, `Midpoint` and `Slope` ; see Section 2.2.2 for the meaning of the parameters. Points correspond to estimates, while segments are 95% asymptotic confidence intervals based on Fisher’s z-transformation. The dashed red line has been added for reference to indicate zero correlation. 34 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.05.764527; this version posted November 6, 2026. The copyright holder for this](https://doi.org/10.64898/2026.11.05.764527) preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. ## 1058 Appendix A Supplementary Figures Figure 11: Trace plots for several model parameters. There are 4 Markov chains per parameter. 35 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.05.764527; this version posted November 6, 2026. The copyright holder for this](https://doi.org/10.64898/2026.11.05.764527) preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 36 Figure 12: The autocorrelation function for Markov chains of several model parameters. The figure is best viewed sideways. There are 4 Markov chains per parameter. Each is identified by a number from 1 to 4. [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.05.764527; this version posted November 6, 2026. The copyright holder for this](https://doi.org/10.64898/2026.11.05.764527) preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. Figure 13: The residuals versus fitted values plots for an arbitrary selection of 4 genotypes in 4 trials. 37 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.05.764527; this version posted November 6, 2026. The copyright holder for this](https://doi.org/10.64898/2026.11.05.764527) preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. ## 1059 Appendix B Supplementary Listings and Tables Listing 1: Variable and parameter names in the **brms** model. # Env : t r i a l i d e n t i f i e r # geno : genotype i d e n t i f i e r # biomass : scaled biomass # dscaled : scaled time # Asym: carrying capacity of the log − l o g i s t i c function # tmid : midpoint parameter of the log − l o g i s t i c function # slope : ( scaled ) slope parameter of the log − l o g i s t i c function 38 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.05.764527; this version posted November 6, 2026. The copyright holder for this](https://doi.org/10.64898/2026.11.05.764527) preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. Listing 2: Model formula as implemented in **brms** . formula brms _<_ − bf ( biomassn l f ( mubiomass˜ ˜ mubiomass,˜ Asym / (1 + ( dscaled / tmid)ˆ( − slope ) ) ), n l f (Asym ˜ exp ( logAsym ) ), n l f ( tmid ˜ exp ( logtmid ) ), logAsymlogtmidl o g s l o p en l f ( slope˜˜ ˜ 11 1 exp ( l o g s l o p e+++(1(1(1 _||_ _|_ Env)Env)Env)) )+++, (1(1(1 _||_ _|_ geno )geno )geno )+++(1(1(1 _||_ _|_ geno : Env ),geno : Env ),geno : Env ), nl = TRUE ) 39 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.05.764527; this version posted November 6, 2026. The copyright holder for this](https://doi.org/10.64898/2026.11.05.764527) preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. Listing 3: Priors for the **brms** model. prior brms _<_ − p r i o r s t r i n g (” normal (0, 2)”, c l a s s = ”sigma ”) + p r i o r s t r i n g (” normal (2, 2)”, c l a s s = ”b”, nlpar = ”logAsym ”, coef = ” I n t e r c e p t ”) + p r i o r s t r i n g (” normal (2, 2)”, c l a s s = ”b”, nlpar = ” logtmid ”, coef = ” I n t e r c e p t ”) + p r i o r s t r i n g (” normal (2, 2)”, c l a s s = ”b”, nlpar = ” l o g s l o p e ”, coef = ” I n t e r c e p t ”) + p r i o r s t r i n g (” exponential (1)”, c l a s s = ”sd ”, group = ”geno ”, nlpar = ”logAsym ”) + p r i o r s t r i n g (” exponential (1)”, c l a s s = ”sd ”, group = ”geno ”, nlpar = ” logtmid ”) + p r i o r s t r i n g (” exponential (1)”, c l a s s = ”sd ”, group = ”geno ”, nlpar = ” l o g s l o p e ”) + p r i o r s t r i n g (” exponential (1)”, c l a s s = ”sd ”, group = ”Env”, nlpar = ”logAsym ”) + p r i o r s t r i n g (” exponential (1)”, c l a s s = ”sd ”, group = ”Env”, nlpar = ” logtmid ”) + p r i o r s t r i n g (” exponential (1)”, c l a s s = ”sd ”, group = ”Env”, nlpar = ” l o g s l o p e ”) + p r i o r s t r i n g (” exponential (1)”, c l a s s = ”sd ”, group = ”geno : Env”, nlpar = ”logAsym ”) + p r i o r s t r i n g (” exponential (1)”, c l a s s = ”sd ”, group = ”geno : Env”, nlpar = ” logtmid ”) + p r i o r s t r i n g (” exponential (1)”, c l a s s = ”sd ”, group = ”geno : Env”, nlpar = ” l o g s l o p e ”) 40 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.05.764527; this version posted November 6, 2026. The copyright holder for this](https://doi.org/10.64898/2026.11.05.764527) preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. Table 5: Summary of the **brms** fit. ``` Family: gaussian Links: mu = identity; sigma = identity Formula: biomass ~ mubiomass mubiomass ~ Asym/(1 + (dscaled/tmid)^(-slope)) Asym ~ exp(logAsym) tmid ~ exp(logtmid) slope ~ exp(logslope) logAsym ~ 1 + (1 | Env) + (1 | geno) + (1 | geno:Env) logtmid ~ 1 + (1 | Env) + (1 | geno) + (1 | geno:Env) logslope ~ 1 + (1 | Env) + (1 | geno) + (1 | geno:Env) Data: data_small_subsampled (Number of observations: 61213) Draws: 4 chains, each with iter = 2500; warmup = 500; thin = 1; total post-warmup draws = 8000 Group-Level Effects: ~Env (Number of levels: 124) Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS sd(logAsym_Intercept) 0.42 0.03 0.37 0.48 1.01 414 986 sd(logtmid_Intercept) 0.33 0.02 0.29 0.38 1.00 442 742 sd(logslope_Intercept) 0.16 0.01 0.14 0.18 1.00 648 1216 ~geno (Number of levels: 48) Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS sd(logAsym_Intercept) 0.07 0.01 0.05 0.08 1.00 1390 2427 sd(logtmid_Intercept) 0.05 0.01 0.04 0.06 1.00 1365 2506 sd(logslope_Intercept) 0.04 0.00 0.03 0.05 1.00 2029 3698 ~geno:Env (Number of levels: 5952) Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS sd(logAsym_Intercept) 0.05 0.00 0.05 0.05 1.00 2860 5605 sd(logtmid_Intercept) 0.01 0.00 0.01 0.01 1.00 1764 3395 sd(logslope_Intercept) 0.00 0.00 0.00 0.00 1.00 2708 3828 Population-Level Effects: Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS logAsym_Intercept 2.21 0.04 2.14 2.28 1.01 186 419 logtmid_Intercept 2.09 0.03 2.02 2.15 1.03 119 182 logslope_Intercept 1.29 0.01 1.26 1.32 1.01 303 653 Family Specific Parameters: Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS sigma 0.18 0.00 0.17 0.18 1.00 5237 4733 Draws were sampled using sample(hmc). For each parameter, Bulk_ESS and Tail_ESS are effective sample size measures, and Rhat is the potential scale reduction factor on split chains (at convergence, Rhat = 1). ``` 41