where wi is the relative fitness of individual i. The µw denotes the overall fixed effect mean fitness, and CohortFi a fixed effect denoting the year of first breeding of individual i, which is entered to correct for temporal variation in fitness and the possible truncation of lifetime fitness in the last cohorts considered in the analysis. The bivariate element in the random regression is specified by the polynomial . Thus, for order x, the variances in and are both estimated, in addition to all possible covariances between these parameters. Nevertheless, because each individual has only one estimate of lifetime fitness, only the variance in elevation ( ) is estimable and all other variances (including the residual variance) must be constrained to zero. Thus, for a first-order random regression, the mixed model procedure needs to estimate the (co)variances (written in matrix form, omitting the upper diagonal that is symmetric to the lower diagonal).
where denotes the selection differential on the elevation of the clutch size–laying date relationship (Selevation) and the selection differential on the slope of the clutch size–laying date relationship (Sslope). These selection differentials can be expressed as a selection gradient by dividing them with and , respectively, or as selection intensity (standardized selection gradient) by dividing them with the square root of these variances.
Higher order polynomials can be introduced as an extension of this matrix. Thus, bivariate random regression allows estimating directional selection on the properties of a polynomial individual-specific relationship between labile life-history traits. Testing the formal significance of the selection differential can be performed by constraining the appropriate covariance to zero and performing a likelihood ratio test between the unconstrained and constrained models. The procedure does not, however, allow description of nonlinear selection. In order to evaluate the potential for nonlinear selection on these properties, we suggest investigation of the map of relative fitness on the BLUP values for all reaction-norm properties. Under the paradigm of optimal reaction norms, one would expect that variation in slope may be under stabilizing selection, which can hereby be graphically evaluated.
All random regression models were solved using Restricted Maximum Likelihood (REML) in the programe ASReml (VSN International). This programe uses the delta method (see e.g., Lynch and Walsh 1998 ) to estimate the standard errors of functions of variances. The code used in this paper and example file of data is provided in the Supporting Information. In the implementation of the above equations in AsReml, the following need to be observed. (1) In random regression with a linear covariate, the variances of first-order and higher order terms are dependent on the scaling of the covariate (e.g., Pinheiro and Bates 2000 ; Schaeffer 2004 ), which is dependent on the software used to implement the models. AsReml standardizes the covariate such that its minimal value is –1 and the average of the covariate values is zero. Hence, variances in elevation are defined for the average of all the covariate values (which may or may not equal the average covariate when taking the number of observations into account). Random-regression slope is then defined as change in response variable per unit equal to the difference between minimum and mean value of the covariate. (2) The within-individual variance and variance in slope(s) of lifetime fitness cannot be constrained to exactly zero and must, in practice, be constrained to a very small value. See the Supporting Information for additional details.