Microevolutionary system identification and climate response predictions by use of BLUP prediction error method

Animals and other organisms in wild populations may adjust to climate change by means of plasticity and evolution, and it is an important task to find the contributions from each of these effects. Attempts to solve this disentanglement problem by use of best linear unbiased prediction (BLUP) and restricted maximum likelihood (REML) methods, as borrowed from the field of domestic breeding, have been criticized because of errors in the variances of the predicted random effects. A primary purpose of this article is to show how the problem can be solved by use of BLUP in a prediction error method (PEM), borrowed from the well-established engineering system identification discipline. The PEM approach is first to collect environmental input data u t and mean phenotypic output data y t , as well as individual phenotypic and fitness data, for consecutive generations from t = 1 to T . A reaction norm model of the evolutionary system is then used to find predictions (cid:98) y t , and the parameters in this model, together with environmental reference values and initial state variables, are finally tuned such that (cid:80) Tt =1 (cid:16) y t − (cid:98) y t (cid:17) 2 is minimized. The main contribution is the use a dynamical BLUP model in a BLUP/PEM method for parameter estimation and mean reaction norm trait predictions. The model is dynamical in the sense that the incidence matrix in an underlying linear mixed model, as well as the corresponding residual covariance matrix, are functions of time. For comparisons, a selection gradient prediction model as presented in Ergon (2022a,b) is also used in a GRAD/PEM method. The advantages of the BLUP/PEM method are that it can utilize genetic relationship information, and that it produces better estimates of environmental reference values. The treatment is limited to multiple-input single-output (MISO) systems. Generations are assumed to be non-overlapping. Simulation examples show that BLUP/PEM may find good estimates of environmental reference values and initial state variables, as well as good mean reaction norm trait predictions. Details for use of additional fixed effects, as well as appropriate methods for model validation remain to be worked out.


Introduction
Individuals in a population of wild organisms may adjust to climate change by means of phenotypic plastic-ity, i.e., by changes in individual behavioral, morphological and physiological traits in response to changes in the environment, without natural selection being involved.A population of wild organisms may also adapt to climate change by evolution, i.e., by natural selection such that the proportion of individuals with favourable traits increases from generation to generation.The combined effect of plasticity and evolution can, for example, be seen in that many bird species have advanced their breeding period by several weeks in the last decades (see, e.g., Bowers et al., 2016).It is well documented that individual phenotypic plasticity plays an important role in such acclimations (Merilä and Hendry, 2014), but in many cases it is unclear to which extent evolution plays a role, and this is further complicated by the fact that also plasticity may evolve.Methods for disentanglement of plasticity and evolutionary effects on acclimation to climate change is therefore an active research area (see, e.g., Scheiner et al., 2020).A main purpose of this article is to show that system identification methods from the field of engineering control can play an important role in such disentanglement efforts.
The combined effect of plasticity and evolution is often modeled by use of reaction norms, for example the two-trait intercept-slope model as used in Lande (2009), where y i,t is the individual phenotypic response on an environmental cue u t − u ref at generation t (where Lande assumed u ref = 0 and η i,t = 0), and where u ref is a reference environment as discussed below.Note that I in the following will treat the reaction norm parameters a i,t + v i,t and b i,t + η i,t as quantitative traits in their own right, and where necessary I will thus distinguish between phenotypic traits and reaction norm traits.In Equation (1a), a i,t and b i,t are the additive genetic components of the reaction norm traits, while v i,t and η i,t are independent and identically distributed (iid) zero mean normal non-additive effects (genetic or environmental white noise).From Equation (1a) follows the mean trait equation where the mean traits a t and b t may evolve as results of natural selection (Lande, 2009).The plasticity-evolution disentanglement problem for a population system given by Equations (1a,b) is illustrated in Figure 1, where an individual Gaussian fitness function is introduced.Fitness is here a measure of reproductive success (Ergon, 2019(Ergon, , 2022a)), for example number of offspring.Initially, it is assumed that the population is fully adapted to the temperature u ref = 10 • C, in the sense that the expected geometric mean fitness has a global maximum at this temperature, which implies that the variance of y i,t has a global minimum (Lande, 2009).Panel b) shows 25 realizations of individual reaction norms, as well as the mean reaction norm y t , in a stationary stochastic environment with the expected temperature 10 • C, and thus the expected cue value E [u t − u ref ] = 0.As shown in panel a), this means that E [y t ] = 0 is equal to the expected phenotypic value that maximizes fitness, such that the variance in individual fitness values is at a minimum.If the temperature then suddenly changes to 15 • C, the expected mean phenotypic value will immediately change to E [y t ] = 1.2, and this will be a purely plastic response.If at the same time the peak of the fitness function is assumed to move from 0 to 2, as shown in panel c), it follows that the mean fitness can be improved by evolution of both a t and b t , such that E [y t ] finally reaches the value 2. Note that the mean fitness will then be lower than the initial mean fitness, because the variance in individual fitness values will be larger.Concomitant with a settlement on E [y t ] = 2, the value of E b t will return to the initial optimal value in a stationary stochastic environment (Lande, 2009).In the end, 100% of the change in E [y t ] will thus be due to evolution of a t , while 0% is due to evolution of b t , but evolution of b t may still play a vital role in the period before settlement in the new environment.Evolution of b t may in fact be crucial for survival of the species.See also Lande (2009) for simulations of the mean reaction norm trait responses in a case as described.In the example as illustrated in Figure 1, there is the possibility that the temperature u ref which the population is fully adapted to also is an evolvable trait, in addition to the mean intercept and mean plasticity slope.Such cases are discussed in Ergon and Ergon (2017) and Ergon (2018).
The first essential question in this article is how the evolution of the mean reaction norm traits a t and b t in a population can be found from recorded values of u t , y i,t , and individual fitness, and then also in more realistic cases with a gradual change in a multivariate environment.The second essential question is how a model for this purpose can be identified in a way such that genetic relationships between the individuals in the population are taken into account.As discussed below, also the reference environment and initial mean reaction norm trait values may be identified from time series data.
The adaptive contribution from evolution can be predicted by use of a linear mixed model with fixed and random effects, combined with BLUP (best linear unbiased predictions) and REML (restricted maximum likelihood) methods, as extensively used in domestic breeding (Henderson, 1950;Ch. 26 and 27 Lynch and Walsh, 1998).BLUP/REML applications on wild animal populations are for example discussed in Kruuk (2004) and Nussey et al. (2007), while Arnold et al. (2019) discusses plant applications.Hadfield et al. (2010) criticized this approach on the ground that BLUP underestimates the variances of the random effects in the linear mixed model.However, as shown in Ergon (2022c), this underestimation is precisely what is needed, in order to obtain the correct incremental changes for prediction of phenotypic traits of non-plastic organisms, and as shown in Ergon (2023), this is the case also for reaction norm traits of plastic organisms.
With known u t and y t time series data, the traditional BLUP/REML breeding problem is turned into what in the engineering control community is called a system identification problem.As will be shown in Section 2, the use of input data u t in a linear mixed model leads to a time varying coincidence matrix Z Z Z t , and REML at a single generation cannot therefore be used for parameter estimation, as it is otherwise done for BLUP parameter estimation.An alternative is therefore to apply REML on data from several generations, or to use a prediction error method (PEM) that utilizes data from all generations in a study.System identification can be described as the art and science of building mathematical models of dynamical systems from observed input-output signals, and it has a wide range of applications in various fields of physics and engineering (Ljung, 2010;Ljung et al., 2011).PEM is one of the important methods in system identification (Ch. 7-9, Ljung, 1999;Ch. 7, Stoica and Söderström, 1989), and a primary purpose of this article is to show that BLUP/PEM is a feasible solution to the problem with a time varying design matrix in a linear mixed model.Note that both REML and PEM make use of various versions of the numerical Newton algorithm, and that there are clear theoretical links between ML and PEM (Ljung, 2010;Ch. 7, Stoica and Söderström, 1989).
I will thus show how the evolutionary system identification problem can be solved by use of BLUP models in PEM.The treatment will be limited to the relatively simple cases where a population evolves in response to changes in a physical multivariate environment (microevolution).The more realistic cases with evolutionary interaction between several species (macroevolution) may be far more involved.I will discuss evolutionary systems with multiple environmental input signals and a single mean phenotypic output signal (MISO systems), i.e., extensions of the model in Equations (1a,b).I will thus assume an individual multivariate interceptslope reaction norm model where y i,t is the individual phenotypic response on en- define the reference environment.Note that such reference values are needed for all interval scaled input variables, as for example temperature in • C, and that they must be properly defined (Ergon, 2022a).Ratio scaled input variables have a well-defined zero-point, but it follows from the discussion in Ergon (2022a) that these zero-points are not necessarily the proper reference values.With multiple input variables, the reference environment is defined by a column vector u u u ref .
In Equation (2a), a i,t , b 1,i,t , b 2,i,t , . . ., b q,i,t are the additive genetic components of normally distributed individual reaction norm intercept and plasticity slopes, respectively, and it is thus the mean values a t , b 1,t , b 2,t , . . ., b q,t that are evolvable.The variables a i,t , b 1,i,t , b 2,i,t , . . ., b q,i,t may be correlated, as determined by an additive genetic covariance matrix G G G. The terms v i,t , η 1,i,t , η 2,i,t , . . ., η q,i,t are iid zero mean normal non-additive effects (genetic and environmental white noise).As done in Lande (2009) and Ergon and Ergon (2017), the variables a i,t +v i,t , b 1,i,t +η 1,i,t , b 2,i,t +η 2,i,t , . . ., b q,i,t + η q,i,t are treated as quantitative traits in their own right.These traits have a covariance matrix P P P = G G G + H H H, where H H H is a diagonal matrix with variance values σ 2 v , σ 2 η1 , σ 2 η2 , . . ., σ 2 ηq along the main diagonal.The time unit is generation, and generations are assumed to be non-overlapping.
From Equation (2a) follows the equation for the mean phenotypic value, Here, the mean reaction norm traits may be collected . The theory will be developed under the assumption of constant matrices G G G and P P P , but a case with time-varying matrices will be tested in simulations.
Although the main purpose is to show how a BLUP model can be incorporated into the PEM approach, I will also compare with results found by a selection gradient (GRAD) method as introduced in Ergon (2022a), and as further developed in Ergon (2022b).This method is based on the multivariate breeder's equation (Lande, 1979).Principal differences between the two approaches will be further discussed below.
Independent of the details of the PEM identification method, there are some general system identification problems that must be considered.First, the initial state of the system affects the dynamics, and this state must therefore be known or identified (Liu, 2013).A complete identification of an unknown dynamical reaction norm system will thus require that also the initial mean reaction norm traits are estimated.
Second, PEM uses a mathematical model with the same input signal as the true system, such that the system can be identified by parameter optimization (model tuning) until the model output in some sense is as close as possible to the true output.For reaction norm models this implies that also the environmental reference vector u u u ref must be known or estimated, and that mean centering of for example a temperature trend is not a good idea (Ergon, 2022a).Such mean centering appears to be common in applications of linear mixed models in evolutionary contexts, but that assumes stationary stochastic environments (e.g., Thomson et al., 2018).In practice it may be difficult to find the reference point in the environment space, and it is even more difficult because it may evolve (Er-gon and Ergon, 2017;Ergon, 2018).I will, however, assume a constant reference environment.
Third, the aspect of system identifiability is essential, and an identifiability analysis of the GRAD and BLUP microevolutionary models is therefore a part of Section 2, including identifiability of the reference environment and initial mean reaction norm trait values.Here, note that traditional REML as used in domestic breeding (Ch. 27, Lynch and Walsh, 1998) would have to be modified for estimation of the reference environment.As verified by simulations, such an estimation is possible in the BLUP/PEM approach.
As we will see, microevolutionary PEM relies on measurements of the environmental input variables, as well as on individual phenotypic values and fitness, and modeling errors and noise corrupted measurements can make the estimation of especially u u u ref difficult.As explained in Section 2, and shown in simulations, this is a serious problem in the GRAD/PEM approach, while the BLUP/PEM approach is much less sensitive.A main question is still how errors in u u u ref affect the mean reaction norm trait predictions.
An additional aspect is that although dynamical systems to some extent can be identified from output-only measurements, it will in practice often be impossible to find good parameter estimates unless the system is persistently excited by the input signals (Ch. 13, Ljung, 1999).This is hardly a problem when we as in the simulations in Section 3 assume yearly temperatures in for example wild bird habitats as part of the environmental input signals, and persistent excitation of sufficient order will therefore be assumed.
The PEM principle is simple, as illustrated in Figure 2.For readers who are not accustomed to block diagrams, it should be pointed out that the block with an adjustable dynamical system model (the tuning model) is just an illustration of either the GRAD or dynamical BLUP prediction equations, as developed in Section 2. These equations make use of environmental input data u u u t , and phenotypic output data y y y t (a vector of individual values) and y t (the population mean value) from the microevolutionary system, which are collected for a number T of consecutive generations.In addition, a vector w w w t of individual relative fitness values must be collected for each generation, where the individual fitness is a measure of reproductive success.The relative fitness is simply the individual fitness divided by the mean fitness in the population under study.
Note that the individual fitness values are used as inputs to the tuning model, and that the tuning model therefore does not include a fitness function.In simulations, however, an individual fitness function is needed in the evolutionary model that generates true data.As illustrated in Figure 1., this function must define the phenotypic value that gives maximum fitness as function of the environment.
The inputs u u u t , y y y t and w w w t to the tuning model (the GRAD or BLUP prediction equations) result in the predicted mean phenotype y t , such that the prediction error y t − y t can be computed.In BLUP/PEM, also a known additive genetic relationship matrix A A A t is necessary (Ch. 26, Lynch and Walsh, 1998).The model parameters, including u u u ref and the initial mean reaction norm traits, are finally adjusted such that is minimized.The art part of this solution is to formulate a tuning model that does not have unnecessarily many parameters to estimate, while it at the same time is rich enough to capture the essential properties of the microevolutionary system.
T .Here, u u u t and y t are the known environmental input vector and the known mean phenotypic value at time t, respectively.A A A t is the known additive genetic relationship matrix, while y y y t and w w w t are known vectors of individual phenotypic and relative fitness values, respectively.The G G G and P P P matrices include the system parameters, while x x x init and u u u ref are vectors defining the initial mean reaction norm trait values and the reference environment, respectively.All of these model parameters are tuned until The microevolutionary prediction error method as illustrated in Figure 2 will, as already mentioned, be investigated with use of two types of tuning models, where predictions y t = a t + (u u u t − u u u ref ) T b b b t are found in different ways.Primarily, y t is found from a dynamical BLUP model, which makes use of genetic relationship information.For comparisons, y t is also found by use of a selection gradient (GRAD) model, as developed for a SISO system in Ergon (2022a), which ignores information on genetic relationships in the population.
After tuning of the prediction model in Figure 2, i.e., when y t ≈ y t for all generations from t = 1 to T , predictions of the mean reaction norm traits can be recovered from the model.The quality of these predictions with random measurement errors in the u u u t and y y y t data will here be studied by means of simulations (while w w w t will be the relative number of offspring, which will be assumed to be correctly counted).Prediction errors caused by time-varying parameter values and modeling errors will also be studied.
Section 2 defines the covariance structure for Equation (2a), and gives the theoretical background for the alternative GRAD and BLUP optimization algorithms used for tuning according to A final summary with discussion is given in Section 4. Proof of a theorem is given in Appendix A, while Appendix B shows GRAD results with population size N = 10, 000.MATLAB simulation code is finally given in Appendix C.

Covariance structures
For the GRAD and BLUP modeling it is necessary to define the covariance matrices G G G and P P P = G G G + H H H to be estimated according to Figure 2.For the BLUP modeling it is in addition necessary to introduce Kronecker product covariance matrices G G G t and R R R t .
First, the reaction norm traits in Equation ( 2a) have an (1 + q) × (1 + q) additive genetic covariance matrix given by the expectation The reaction norm traits also have the phenotypic covariance matrix Second, a population of N genetically related individuals has the N (1 + q) × N (1 + q) Kronecker product covariance matrix G G G t = G G G A A A t , where A A A t is the symmetric N × N additive genetic relationship matrix (Ch. 26, Lynch and Walsh, 1998;Cantet et al., 2022;Mathew et al., 2018), and where is the Kronecker product operator.The Kronecker product simply means that each of the elements in G G G are multiplied by A A A t .If the individuals are genetically unrelated, the additive genetic relationship matrix is a unity matrix, A A A t = I I I N .In BLUP modeling, a population of N individuals also has a residual Kronecker product covariance matrix R R R t = R R R t I I I N , where R R R t is the covariance matrix of the residuals in a linear mixed model.

The GRAD prediction model
The GRAD prediction model is based on the multivariate breeder's equation (Lande, 1979), which means that it implicitly assumes that the additive genetic relationship matrix is a unity matrix (Ergon, 2023).Applied on the reaction norm traits in Equation ( 2a), the multivariate breeder's equation for the incremental changes in mean traits from generation to generation becomes where w i,t i th relative individual fitness, i.e., w i,t = W i,t /W t .Note, however, that Equation (3) cannot be used in the tuning model in Figure 2, because the individual traits a i,t , etc., are not available.
In Ergon (2022a), this problem was solved by use of a linear transformation of the vector a where ∆u where Here, with Note that Equation ( 5a) is a dynamical model in the sense that P yy,t varies with time, and that we can set G aa to any fixed value and scale the rest of the parameters accordingly.Also note that Equations ( 3) and (5a) give identical results only asymptotically, i.e., for N → ∞.In order for Equation (5a) to be compatible with the Price equation (Price et al., 1970;Ergon, 2019), the covariance expression should be computed as With given initial mean trait values, Equation (5a) gives a t , b 1,t , etc., and thus also

Dynamical BLUP prediction model
In order to handle cases where the additive genetic relationship matrix A A A t departs from a unity matrix, we may use a linear mixed model and BLUP solutions.As a starting point we may note that the reaction norm model in Equations (2a,b), with use of , and a i,t − a t = a ′ i,t etc., can be combined into We thus use the individual variations around the mean values a t , b 1,i,t , . . ., b q,t as individual random effects.This leads to the linear mixed model for a population of N individuals, where ] and E η η η j,t in this mixed model are zero by definition (Ch. 26, Lynch and Walsh, 1998;Robinson, 1991).From Equation (6b) follows the residual covariance matrix R R R t = E e e e t e e e T t = σ In order to predict y t , a a a ′ t , and b for any given parent generation, we may solve the BLUP equation (Ch. 26, Lynch and Walsh, 1998;Robinson, 1991;Henderson, 1950) Note that this is a dynamical model in the sense that Z Z Z t , G G G t and r t vary with time.
T , and since G G G , it follows from Equation (6e) that the BLUP equation for any given parent generation is where Here, the predictions are unbiased in the sense that E y t = E [y t ] and E a a a ′ t = 0 0 0, E b b b ′ 1,t = 0 0 0, etc. (Ch.26, Lynch and Walsh, 1998;Robinson, 1991).From Equation (6f) follows the following theorem: Theorem 1 In Equation (6f), G aa may be set to any fixed value, if at the same time the rest of the parameters in G G G as well as in r t are given normalized values in relation to the value of G aa .
Proof : The term r t G G G −1 in Equation (6f) remains constant if all elements in G G G and r t are multiplied by the same factor.
Next, we must find how fitness affects the mean reaction norm trait predictions for the offspring of any given parent generation.As shown in Ergon (2023), and as verified by Theorem 2 below, the incremental changes in mean trait values can be found from the predicted random effects by use of Robertson's secondary theorem of natural selection (Ch.6, Walsh and Lynch, 2018;Robertson, 1966) where a ′ i,t , b ′ 1,i,t , . . ., b ′ q,i,t are found from Equation (6f).A comparison with the multivariate breeder's Equation (3) indicates that Equation ( 7) assumes that G G G = P P P , which in general is not correct.However, as discussed in Ergon ( 2023 7), without direct use of u u u ref .

Comparison of prediction methods for
A A A t = I I I N Under the assumption that A A A t = I I I N , it follows from Theorem 2 in Ergon ( 2023) that the BLUP model described above leads to the mean trait prediction equation where . A comparison of the GRAD and BLUP prediction models leads to the following theorem, which follows from Theorem 3 in Ergon ( 2023): Theorem 2 With A A A t = I I I N , the GRAD and BLUP models result in exactly the same mean reaction norm trait predictions, i.e., Equation (5a) and Equation (8) give identical results, independent of population size.

Identifiability and prediction error issues
Ideally, optimization according to Figure 2 gives y t ≡ y t for t = 1 to T , and estimated entities u u u ref ≡ u u u ref , , and P P P ≡ P P P .That would require an input u u u t that is persistently exciting of sufficient order (Ch.13, Ljung, 1999), a perfect tuning model with perfect measurements of the input and output signals, and that T → ∞.In order for the system to be parameter identifiable (Ch.6, Stoica and Söderström, 1989), it should in theory be possible to obtain y t ≡ y t for t = 1 to T only with u u u ref ≡ u u u ref , It may be argued that P yy,t according to Equation (5b) includes the term G aa + σ 2 v , which indicates a structural identifiability problem.This is, however, not the case when G aa is given a fixed value in accordance with Theorem 1.
The identifiability property of the GRAD system in Equations ( 2a) and (5a) is given by the following theorem: Theorem 3 The system given by Equation (2a) and the GRAD Equation (5a) is parameter identifiable, provided that G aa ̸ = 0 exists, and that G aa is given a fixed value.
See Appendix A for detailed proof.
It follows from Theorems 2 and 3 that also a BLUP system in the special case with a genetic relationship matrix A A A t = I I I N is parameter identifiable.There seems to be no reason why the identifiability properties of a BLUP system should be affected when A A A t ̸ = I I I N , and this reasoning leads to the following conjecture, which is also supported by simulation results in Section 3: Conjecture 1 Theorem 3 is valid also for the system given in Equation (2a) and the BLUP Equation (6f).
When the reference environment is fixed but incorrect, the situation is different.It follows from Equation (5a) that in order to minimize T t=1 y t − y t 2 , the effects of the error in u u u ref must somehow be compensated by errors in G G G or P P P .The resulting deviations in parameter values may be large, and as shown by simulations in Section 3 this will result in errors in especially the predicted mean intercept value.As discussed in Ergon (2022a), the effects on the predicted mean slope values are smaller.

Simulation Results
I will here assume that true data is generated by a BLUP model, and test the GRAD and BLUP optimization methods in several simulation examples.In Example 1 we will see that the GRAD tuning model has problems regarding estimation of the environmental reference value, and that is the case also if true data is generated by a GRAD model.In Example 2 we will see that the BLUP tuning model gives good estimates of u ref , and good predictions of the mean reaction norm traits, also with measurement errors.With perfect measurements, also the parameter estimates will be good.This example also shows that a fixed but incorrect value of u ref may give large prediction errors, and it verifies that the BLUP and GRAD predictions with an additive genetic relationship matrix A A A t = I I I N are identical (Theorem 2).Example 3 shows that BLUP optimization may give good mean trait predictions also with true data generated by use of time-varying parameters, but then with large errors in the parameter estimates.Finally, Example 4 verifies that the final value of the criterion function T t=1 y t − y t 2 may be used for model selection (see discussion in Section 4).

Toy example, fitness function and input signals
For the purpose of testing the BLUP/PEM optimization method, I assume a true system with reaction norms according to Equations (2a,b), where y i,t is the individual clutch-initiation date for a certain bird species.This is similar to the toy system in Ergon (2022a), except that a second input signal is added, and that I will assume a continuous clutch-initiation time (which may be more realistic).As in Ergon (2022a), the individual fitness is an integer number from 0 to 10 (the number of offspring).
The input signals are the present spring temperature u 1,t = u t , and the spring temperature u 2,t = u t−1 a generation back in time.The delayed effect could for example be food abundance in year t determined by the temperature in year t − 1.This means that the two input signals have the same reference environment, and for simplicity I assume a temperature scale such that u ref = 0.
The individual fitness values are rounded values of where θ t is the phenotypic value that maximizes fitness, while ω 2 = 1000 is the squared width of the Gaussian fitness function.Assume a constant or slowly varying mean µ U,t of a stochastic environment, with added white noise (iid zero mean normal random variations) u n,t , with variance σ 2 Un = 0.5, i.e., u t = µ U,t + u n,t .In a corresponding way assume that θ t = µ Θ,t + θ n,t , where µ Θ,t = −20µ U,t and θ n,t is white noise with variance σ 2 Θn = 200.Also assume that u n,t and θ n,t are correlated with covariance cov (u 1,t , θ t ) = −2.5, as explained from fluctuations in the environment during the development of the individuals in the parent generation (Lande, 2009).Since u t is white noise and u 2,t = u t−1 , we have cov (u 2,t , θ t ) = 0.
Further assume that u t and θ t are noisy ramp functions as shown in Figure 3, starting at t = 10.The choice of a negative trend in θ t , and thus in y t , is inspired by common cases where a positive temperature trend leads to earlier breeding dates for various natural populations (e.g., Bowers et al., 2016).The noisy temperature trend in Figure 3 is similar to the registered yearly mean temperature in Oslo, Norway, from 1960 to 2020 (klimaservicesenter, 2023), using the temperature in 1960 as zero-point.In order to make the simulation examples more realistic, only data from t = 31 to 60 (1991 to 2020) are used for parameter estimation and mean reaction norm trait predictions.This will in most realizations mean that the reference environment is not within the range of input data used in the optimizations.

Population size, measurement errors, initial values, and result presentation
The GRAD and BLUP optimization algorithms are tested in four examples in Subsection 3.4 below, with in all 10 special cases.The population size is in general N = 100, although for comparison a test with N = 10, 000 is included in Appendix B. The algorithms are tested with perfect input and output information, as well as with measurement errors in u 1,t , u 2,t and y i,t .The errors in temperature measurements are implemented by adding normally distributed white noise with standard deviations 0.1 to u 1,t and u 2,t , which is 14% of the standard deviation in the yearly variations.The errors in individual phenotypes are implemented as uniformly distributed random errors in the interval −0.5 to 0.5, i.e., errors up to half a day in clutchinitiation time.The relative individual fitness w w w t , i.e., the relative number of offspring, will be assumed to be correctly counted.In all examples the true value G aa = 5 was used, such that estimates of the other G G G and P P P parameters are found relative to G aa (Theorem 1).The initial values for these parameters were sat to 0. The initial mean trait values were sat to b 1,init = 0 and b 2,init = 1 (because the true value is 0), while a init follows from Equation (2b) with y init set to zero.In cases where u ref is a free variable, the initial value was sat to 1 (because the true value is 0).The GRAD optimization with population size N = 100 takes around 1 second on an HP EliteBook ×360 1030 G3 laptop.With N = 100, the optimization time for the BLUP optimization increased to around 400 seconds, owing to the repeated inversions of the (3N + 1) × (3N + 1) matrix in Equation ( 6f).No attempts were made to speed up this optimization.
Simulation results are presented as plots over predicted mean values y t , a t , b 1,t and b 2,t , as compared to true mean values y t , a t , b 1,t and b 2,t .Estimated parameter values are given in tables as mean values and standard errors, M ean ± SE, based on 10 repeated simulations with different realizations of random inputs (in Case 6 increased to 100).The relative errors in total change of predictions over 30 generations are included in the tables, computed as, for example, ∆ error 30 a t % = 100 ∆ 30 a t − ∆ 30 a t /∆ 30 a t etc., where ∆ 30 a t = a 60 − a 31 and ∆ 30 a t = a 60 − a 31 .The final values ε 2 t,f inal of 60 t=31 ε 2 t are also included, as they may possibly be used for model selection.

Generation of true BLUP data
True mean responses y t , a t , b 1,t and b 2,t are generated by means of the dynamical BLUP model (6f) and the prediction Equation ( 7), the fitness function in Equation (9), and known parameter values in the covariance matrices Here, G ab1 = 0 in the true system, but left as a free parameter in the optimizations.It is assumed that the population is fully adapted to the stationary stochastic environment up to t = 10 (1970), i.e., to the expected spring temperature E [u t ] = 0 (Ergon, 2022a).From this follows that G ab1 = G ab2 = 0 (Lande, 2009)), and that the optimal mean plasticity slope values in a stationary stochastic environment are b 1,init,true = cov (u 1,t , θ t )/σ 2 Un = − 5, and b 2,init,true = cov (u 2,t , θ t )/σ 2 Un = 0 (Lande, 2009;Ergon and Ergon, 2017).This will thus be the initial mean plasticity slope values at generation 1.
New random effects were at each generation found as , and where The different data vectors z z z 0,t were drawn from normal distributions in accordance with the given G G G matrix.For the three values of the factor c in Equation ( 10) used in generations 31 to 60, this reduced the variances of the random effects to approximately 95%, 90% and 75% of the nominal values, respectively.The random effects were also mean centred.

Simulation examples Example 1: BLUP true data with GRAD prediction model
Results for different cases with a BLUP true data model and a GRAD tuning model, are given in Table 1.Note that also with perfect measurements (Case 1), the value of u ref is estimated with a large standard error, resulting in a large SE value for ∆ error 30 a t .This is to some extent understandable, since the true data is generated by a BLUP model, but also with true data generated by a GRAD model we find the poor results u ref = −0.07± 0.29.With use of the known value of u ref (Case 2), the parameters G b1b1 and G b2b2 are estimated with small errors, while the estimated white noise variances, σ 2 v , σ 2 η1 and σ 2 η2 have large errors, owing to the modeling error.The mean reaction norm traits are in this case predicted with SE values up to 5%., i.e., GRAD/PEM finds good predictions at the price of large errors in the estimated noise variances.The known value of u ref combined with measurement errors (Case 3), gives SE values up to 11%.With a factor five increase in the measurement errors, the mean trait SE values were up to 40%, without any indication of bias.
Example 2: BLUP true data with BLUP prediction model Simulation results with parameter estimates and mean trait predictions found by use of the dynamical BLUP algorithm, are shown in Table 2 and Figure 4.It is interesting to note that the environmental reference value is estimated far better than with use of the GRAD algorithm, for reasons as explained in Subsection 2.4.With perfect measurements (Case 4), the result is improved to u ref = 0.001 ± 0.002 , and as a result the parameter estimates without measurement errors are close to perfect.With measurement errors (Case 5), the BLUP results for u ref with N = 100 are in fact just as good as for N = 10, 000 with use of true GRAD data and GRAD optimization (Appendix B).Measurement errors result in large errors in some of the parameters, with consequences as discussed in Section 4, but the mean trait SE values are still limited to 8%.With a factor five increase in the measurement errors, the mean trait SE values increased to around 20%, still without any indication of bias.When in addition to measurement errors the incorrect relationship matrix A A A t = I I I N is used in the BLUP optimization algorithm (Case 6), the mean trait predictions are still quite good, but the price to pay is large errors in the parameter estimates.
Table 2 also shows mean trait results for GRAD predictions based on the BLUP parameter results, i.e., when the genetic relationships are ignored.For Case 4 and 5, this leads to around 10% overestimation of the mean trait changes.For Case 6, i.e., with BLUP parameter estimates based on A A A t = I I I N , the GRAD and BLUP predictions are identical, as follows from Theorem 2, and this is the case for population sizes down to N = 2.
The results in Table 2 assumes estimated values of the reference environment.If we instead use fixed values of u ref , we obtain results as in Table 3.Note the large values of ∆ error 30 a t % for u ref ̸ = u ref , while as discussed in Ergon (2022a) the effects on the predicted mean slope values are smaller.Also note the large prediction errors for u ref = 1, i.e., for approximately the initial value of the available temperature data, which it may be tempting to use as a reference.
Since A A A t ̸ = I I I N according to Equation ( 10) reduces the variances of the additive genetic effects with 5 to 25%, for the different values of the factor c, it must   Table 3: Errors in predicted total relative change of a t , b 1,t and b 2,t over 30 generations as for Case 4 in Table 2, except for fixed reference environments used in the BLUP optimization procedure.
be expected that the responses on the gradual change in environment are slowed down.As a test, this effect on the BLUP results were therefore investigated by recording the true change in intercept value from generation 31 to generation 60, with use of the same set of typical u t and θ t input data in 100 repeated simulations.For A A A t = I I I N in the true model, the results was ∆ 30 a = −3.4± 0.1 days, while A A A t ̸ = I I I N gave ∆ 30 a = −2.9±0.1 days, which shows that the responses are indeed slowed down.

Example 3: BLUP true data with time-varying G G G matrix
Results with a true system as in Example 2, but with a parameter value G aa that increases linearly from G aa,1 = 5 to G aa,60 = 10, are shown in Table 4. Mean trait predictions based on G aa = 5, and otherwise true parameter values, will now underestimate the total change ∆ 30 a t (Case 7).However, predictions based on the constant value G aa = 5 and otherwise optimized parameter values give good mean trait predictions (Case 8), because parameter values are found such that the variation in G aa is very much compensated for.When Case 8 is repeated with measurement errors (Case 9), the mean trait predictions are still fairly good, although the SE values increase markedly.
In this case also the SE values for the parameter estimates increase markedly.In Cases 8 and 9 some of the estimated parameter values are far from the true values.
Example 4: BLUP true data with reduced BLUP prediction model In an additional case, Case 10, the BLUP tuning model was reduced by setting G b2b2 = σ 2 η2 = 0, while it at the same time it had to handle measurement noise as in Case 5.This resulted in u ref = 0.03 ± 0.1, ε 2 t,f inal = 0.32 ± 0.21, ∆ error 30 a t % = 35 ± 17 and ∆ error 30 b 1,t % = 3 ± 19, i.e., markedly increased errors in the predicted changes of mean intercept.Note the large increase in ε 2 t,f inal , as compared to Case 5, and this was the case for all realizations, although to various degree.This indicates that the value of ε 2 t,f inal may possibly be used as a tool for model selection (see discussion in Section 4).

Summary and discussion
In summary, there are two main purposes with this article.One is to show that system identification methods borrowed from the field of engineering control may play an important role in the study of evolutionary biological systems.The second is to show that the prediction error method (PEM) can be combined with the best linear unbiased predictions (BLUP) method, which is extensively used in domestic breeding programs.The theoretical foundation is the fact that a reaction norm model can be formulated as a linear mixed model, where the random effects are the individual variations in intercept and plasticity slope values, while the environmental cues are included in a dynamical incidence matrix Z Z Z t .From this follows the additional observation that restricted maximum likelihood (REML) based on data from a single generation cannot be used for parameter estimation, for the simple reason that each scalar element in the residual covariance matrix in the BLUP model includes several residual variances.In the formulation of the linear mixed model, we assumed a MISO system, with linear reaction norms as functions of environmental cues u 1,t , u 2,t , etc., but there is nothing in the theory that prevents us from use of non-linear reaction norms that are also functions of u 1,t u 2,t , u 2 1,t , u 2 2,t , etc., as discussed in Ergon (2018).A more general model for a MIMO system is developed in Ergon (2023).
Two evolutionary models for use in PEM as illustrated in Figure 2 was presented and developed in Section 2. The selection gradient (GRAD) model assumes that the additive genetic relationship matrix is a unity matrix, A A A t = I I I N , which implies that all individuals in the population are genetically unrelated.The BLUP model, on the other hand, can handle cases with inbreeding in the populations, i.e., with A A A t ̸ = I I I N .In both the GRAD/PEM and BLUP/PEM methods, the additive genetic variance G aa can be set to any value, and the rest of the elements in the G G G and P P P matrices estimated in relation to that value (Theorem 1).When A A A t = I I I N , the GRAD and BLUP optimization algorithms give identical results (Theorem 2).It is shown that the GRAD prediction model is identifiable (Theorem 3), and from Theorem 2 thus follows that the BLUP prediction model is identifiable when A A A t = I I I N .It is conjectured that the BLUP model is identifiable also when A A A t ̸ = I I I N , and this conjecture is supported by simulation results.
Although the theory is limited to MISO systems, an extension of the GRAD algorithm to multiple-input multiple-output (MIMO) systems is demonstrated in Ergon (2022b).Such extensions are possible also for the BLUP algorithm (Ergon, 2023).The applicability of PEM in an evolutionary context is, however, not limited to the algorithms discussed, it should in principle work for any type of parametrized prediction model that includes mean reaction norm traits, or other latent variables.
A main theoretical result is that the environmental reference must be known or estimated, as also pointed out in Ergon (2022a).This was recognized by Scheiner et al. (2020), who in their simulations used "a 2000 generation equilibration period, with the mean phenotypic optimum held constant at 0, to allow the population to reach mutation-selection-drift equilibrium".An important practical result is that the BLUP/PEM optimization method is far better than the GRAD/PEM method when it comes to estimation of the reference environment, at least with a single reference environment.
Independent of optimization method, there is an important principal difference between the traditional BLUP/REML approach, and the proposed GRAD/PEM and BLUP/PEM approaches.In the BLUP/REML method there is no feedback loop via a minimization algorithm as shown in Figure 2, where data from all observed generations are used.Instead, parameters in the G G G and P P P matrices are estimated by use of REML for a given base population with an assumed genetic relationship matrix A A A t = I I I N (Ch.19, Walsh and Lynch, 2018), or alternatively a relationship matrix computed from genomic information, and predictions of mean trait responses are then found forward in time in an open-loop fashion.PEM is a closed-loop method, where parameters are estimated, and mean trait responses predicted by use of information at all generations.It is, however, necessary to know the additive genetic relationship matrix A A A t for all generations, which requires some assumptions regarding the first generation in a study.The assumption of a perfectly known pedigree may also then be unrealistic, and kinship estimation based on molecular markers might be a better alternative (Goudet et al., 2018).It is beyond the scope of this article to discuss parameter estimation by application of REML on Equation (6e).It is however obvious that the various components of r t in Equation (6c) cannot be found by REML at a single generation, for the simple reason that only r t appears in Equation (6e).
It should be noted that the primary goal of BLUP/PEM is not to find correct G G G and P P P parameter values as such, but to minimize errors in predictions of the mean phenotypic values y t = a t + (u u u t − u u u ref ) T b b b t , which requires that the errors in the predicted mean reaction norm trait values a t and b b b t are minimized.As further discussed below, this raises questions regarding model validation.With a correct model, perfect measurements, and long enough data series, however, PEM also leads to near perfect parameter estimates (Table 2, Case 4).Hadfield et al. (2010) and other authors are critical to the use of BLUP for natural populations, for the main reason that the estimated variances of the individual random variables have large errors.As shown in Ergon (2022cErgon ( , 2023)), and as utilized in Equation ( 7) and verified by Theorem 2, these errors are just what is needed in order to find the correct mean trait predictions by means of Robertson's secondary theorem of natural selection.Hadfield et al. (2010) also suggest that a good test of evolutionary change would be to simply compare the mean trait values of the first and final generations, as done in this paper.
The GRAD/PEM and BLUP/PEM optimization algorithms have been tested by simulations of an intercept-slope reaction norm toy system in Section 3, and these tests illustrate the feasibility of PEM in the microevolutionary context.The model includes two environmental inputs, first a noisy ramp function in the present year's temperature over 60 years (1961 -2020), similar to the recorded temperature trend in Oslo, Norway, and second the last year's temperature.Corresponding correlated changes in the individual phenotypic value that maximizes fitness are used in the individual fitness function (Figure 3).Only data from the last 30 generations (1991 -2020) are used in the optimizations, i.e., the reference temperature before 1960 is not included in the data.Although the optimal clutchinitiation time is forwarded by 30 days over 30 years (Figure 3), the true phenotypic changes is only 10 to 15 days (Figure 4).The phenotypic response is thus lagging behind the optimal response, which is typical for ramp responses of dynamical systems with time constants.The change in mean intercept over 30 generations is typically −3 days.
The BLUP results in Figure 4, panels a), b), and c), are based on random effects that are drawn from distributions that are affected by the matrix square root M √ A A A t , with the structure of A A A t given in Equation ( 10), which reduces the variances of the random effects up to 75% of the nominal values.As a result of this, the phenotypic response has a tendency of lagging further behind the optimal response.
It should be noted that the simulated MISO system with two inputs has a single and common reference environment u ref , and that there is a major difference between the two algorithms when it comes to estimation of this reference environment.The GRAD optimization algorithm gives large errors in u ref , and correspondingly large errors in the mean trait predictions (Table 1, Case 1), while the BLUP algorithm gives good estimates of u ref (Subsection 2.4 and Table 2).When the GRAD algorithm is applied on data from small populations it may therefore be best to set u ref to a fixed value, based on a judgement on which environment the population is adapted to in the time period before environmental changes.This may also be the case for MISO and MIMO systems with a multivariate reference environment, although such cases need further invetigation.As shown in Table 3, large errors in u ref give large errors also in BLUP predictions, and in more realistic cases the estimated reference environment values must therefore be used with care.An alternative choice would be to use the initial environmental value in a recorded time series, but as shown in Table 3 that may be a bad idea.
Table 1, Case 2, and Table 2, Case 5, show results with measurement errors in the environmental and individual phenotypic data, and in these cases some of the parameters are estimated with large errors.The changes in mean trait values over 30 generations are still predicted without bias, and with standard errors less than 8%.
Table 2, Case 6, shows results with measurement errors and A A A t = I I I N , i.e., with a large error in the additive genetic relationship matrix, and the changes in mean trait values over 30 generations were still predicted with mean errors less than 3%, and standard errors less than 10%.The reason for these good predictions, in spite of the error in A A A t , is that minimization of 30 t=1 y t − y t 2 leads to parameter values that partly compensate for the error in A A A t , and some of these variables may apparently be far from the true values.The results in Case 6 also verify Theorem 2.
As shown in Table 2, Case 4, the BLUP algorithm results in small errors in u ref , and correspondingly good mean trait predictions.As these predictions are not very sensitive to measurement errors (Cases 5 and 6) or errors in A A A t (Case 6), the BLUP/PEM method appears to be a feasible solution for the microevolutionplasticity disentanglement problem introduced in Section 1, at least in cases with a single reference environment.The BLUP method is also more flexible than the GRAD method, in that additional fixed effects can be included in the model (Ch. 19, Walsh and Lynch, 2018).
The large errors in the estimated parameter values caused by measurement noise (Cases 2, 5 and 6) limits the usefulness of the parameter estimates for other purposes, as for example for finding heritabilities of reaction norm traits (Ch. 19, Walsh and Lynch, 2018).A similar situation occurs when the true parameter values vary over time (Table 4, Cases 8 and 9).
Tables 1, 2 and 4 include mean values and standard errors of the final value of the criterion function.Tests with a reduced model (Example 4, Case 10) show that ε 2 t,f inal increases for all realizations, although to various degree.This indicates that the value of ε 2 t,f inal may be used as a tool for model selection, but here a more detailed study is necessary.
The fact that predicted mean trait values a t and b b b t may be quite good also with large errors in the estimated G G G and P P P parameter values, as discussed above and demonstrated in the simulations, raises difficult questions regarding model validation.In system identification in general, identified models are preferably validated by use of data that have not been used for the identification (Ljung, 2010), and in an industrial context it is often possible to obtain such validation data.In the evolutionary context, however, the data series are normally obtained from cumbersome field observations over a small number of generations for a specific population of for example squirrels (Boutin and Lane, 2014).Independent data over the same number of generations for the same species may therefore be impossible to obtain.A possible validation approach could be to use a part of the data for identification and the rest for validation, but as seen in the simulations the dynamics in the responses are quite different in different time periods.A remaining possibility is some form of cross validation, but details of such an approach remain to be worked out.
Details regarding additional fixed effects in the dynamical BLUP model, various sources of drift, etc., also remain to be worked out.Finally, and most importantly, tests on real world data must also be left for further work.

A. Proof of Theorem 3
In order to prove that the GRAD system given by Equations (2a) and (5a) is parameter identifiable (Ch.6, Stoica and Söderström, 1989) when N → ∞.This implies that the system is parameter identifiable if and only if an error in one of the elements in these vectors and matrices, cannot be compensated by errors in other component, in a way that is valid for all generations from t to T .
From Equation (5a) follows that ∆a t and ∆b b b t are determined by input-output and fitness information at time t.Since u u u ′T t+1 is unknown at time t, it thus follows from Equation (4), ∆y t = ∆a t +u u u ′T t+1 ∆b b b t +∆u u u ′T t b b b t , that an error in ∆a t cannot be compensated by an error in ∆b b b t , or vice versa.It thus suffices to prove that only perfect parameter values can give ∆ a t ≡ ∆a t according to Equation (5a).
From Equation (5a) follows that ∆ a t /cov (w i,t , y i,t ) = G aa + G G G ab u u u ′ t / P yy,t = N t / P yy,t , where cov (w i,t , y i,t ) is known.After setting G aa = G aa , and the elements in the 1 × q matrix G G G ab and the q × q matrix G G G bb to G abj and G kj , respectively, we find

Figure 1 :
Figure 1: Illustration of phenotypic and fitness changes as results of a sudden change in mean environment from 10 • C to 15 • C. Panels a) and c) show the individual fitness functions at 10 • C and 15 • C, respectively.Panel b) shows individual and mean reaction norms in a stationary stochastic environment with the expected temperature E [u t ] = u ref = 10 • C (dotted lines and solid line, respectively), from which follows that the immediate plastic response on a mean temperature change from 10 • C to 15 • C, is a change in expected mean phenotypic value from 0 to 1.2 (circles).

Figure
Figure 2: Block diagram of microevolutionary GRAD/PEM and BLUP/PEM for a MISO system, with dynamical tuning model based on a reaction norm model with mean traits vector x x x t = a t b 1,t • • • b q,t T .Here, Figure 2. A theorem shows that the element G aa = E (a i,t − a t ) 2 in the additive genetic covariance matrix G G G may be set to any value, such that the rest of the parameters in G G G and P P P will be estimated in relation to that value.A second theorem shows that a genetic relationship matrix A A A t = I I I N results in identical GRAD and BLUP mean reaction norm trait predictions.A third theorem shows that the GRAD prediction model, as well as a BLUP prediction model with A A A t = I I I N are identifiable, and as supported by simulations it is conjectured that also a BLUP prediction model with A A A t ̸ = I I I N is identifiable.Simulations in Section 3 show the feasibility of GRAD/PEM and BLUP/PEM in the microevolutionary context.The environmental inputs are noisy ramp functions, similar to temperature trends caused by climate change.True mean phenotypic responses y t , and true mean trait values a t and b 1,t , b 2,t , . . ., b q,t are generated by use of the BLUP prediction equations, with known parameter values (including initial mean trait values, reference environment, and additive genetic relationship matrix).Individual additive genetic and environmental effects are at each generation drawn at random from populations with the given covariance matrices G G G and P P P , modified by the additive genetic relationship matrix A A A t .With a known model structure, perfect measurements of the input environment, individual fitness, and individual phenotypic values, will result in close to perfect BLUP parameter estimates and predictions y t .This makes it possible to study the effects of modeling and measurement errors, and of time-varying parameter values.It also makes it possible to study how ignorance of genetic relationship information affects the parameter estimates and mean reaction norm trait predictions.The population size in the main simulations is N = 100.
) T b b b t according to Equation (2b), for use in the tuning model in Figure 2. Initial values of the elements in b b b t for t = 1 must be found as part of the identification, while a 1 follows from Equation (2b) by setting y 1 = 0. See Section 3 for a simulation example, where GRAD estimates of u u u ref , b b b 1 , and the parameters in G G G and r t , are found by use of the function fmincon in MATLAB.This function finds a constrained minimum of a scalar function of several variables starting at initial estimates, in this case the minimum of T t=1 ε 2 t = T t=1 y t − y t 2 .
), and verified by Theorem 2 below, this error is compensated by errors in the predicted random effects a a a ′ t , b b b ′ 1,t , . . ., b b b ′ q,t .The errors in the predicted random effects are thus just what is needed in order to find the correct predicted incremental changes according to Equation (7).With given initial mean trait values, Equation (7) gives a t , b 1,t , . . ., b q,t , and thus alsoy t = a t +(u u u t − u u u ref) T b b b t according to Equation (2b), for use in the tuning model according to Figure 2. Initial values of the elements in b b b t for t = 1 must also here be found as part of the identification, while a 1 follows from Equation (2b) by setting y 1 = 0. See Section 3 for a simulation example, where BLUP estimates of u u u ref , b b b 1 , and the parameters in G and r t are found by use of the function fmincon in MATLAB.2.4.Differences regarding estimation of u u u ref Simulations show that estimation of u u u ref by means of the GRAD optimization method is difficult, while estimation by means of the BLUP optimization method is easier, at least with a single reference value.The fundamental difference is that the GRAD model makes use of u u u ref both in the computation of ∆a t and ∆b b b t by use of Equation (5a), and in the computation of y t = a t + u u u ′T t b b b t , while the BLUP model gives ∆a t and ∆b b b t by use of Equation (
the function fmincon in MATLAB.

Figure 4 :
Figure 4: Simulation results for BLUP system.Panels a), b) and c) show typical results for Case 4 in Table 2, with A A A t as given in Equation (10) (although for population size N = 100).Panels d), e) and f) show typical results for Case 6, with A A A = I I I N used in the BLUP optimization procedure.True y t values are shown by solid blue lines, while final predictions y t are shown by blue dots.True a t , b 1,t and b 2,t responses are shown by green lines, while BLUP predictions a t , b 1,t and b 2,t are shown by black dots.The corresponding GRAD predictions using BLUP parameters are shown by dotted magenta lines.

Table 1 :
Estimation and prediction results with BLUP true data model and GRAD tuning model.Results are based on 10 simulations with different realizations of all random input variables.The population size was N = 100.Case 1: Perfect values of u 1,t , u 2,t and y i,t , and u ref as free variable.Case 2: Perfect values of u 1,t , u 2,t and y i,t , and u ref = u ref = 0. Case 3: Measurement errors in u 1,t , u 2,t and y i,t , and u ref = u ref = 0.

Table 2 :
Estimation and prediction results by use of the BLUP optimization algorithm.Results are based on 10 simulations with different realizations of all random input variables.The population size was N = 100, and u ref is in all cases a free variable.Case 4: Perfect values of u 1,t , u 2,t and y i,t , with A A A ̸ = I I I N used in the tuning model.Case 5: As Case 4, but with measurement errors in u 1,t , u 2,t and y i,t .Case 6: As Case 5, but with A A A = I I I N used in the BLUP tuning model, and number of realizations increased to 100.

Table 4 :
Estimation and prediction results as in Table2, but with G aa increasing linearly from 5 to 10 over 60 generations.Case 7: Results found by use of G aa = 5 and otherwise true parameter values.Case 8: Results found by minimization of N t =