Using the bivariate mixed model [5], Blangero [6] derived the following expression for the genetic variance in phenotype response to two different environments:

where

and are the additive genetic variances of the trait in Environments 1 and 2, in this case two ages, is the additive genetic variance of the phenotypic difference between environments (i.e., trait response), and ρ_{G} is the genetic correlation between environments, in this case ages. The term environment, denoted by E, is used in accordance with standard statistical genetic theory [6]. There is no G × E interaction (i.e., = 0) if = = and ρ_{G} = 1 [6]. To model = 0, and ρ_{G} are parameterized as functions of age, which is the environment of interest in our analyses:

where

is parameterized as an exponential function [7] and ρ_{G} as exponential decay [8], and α, β, and λ are parameters to be estimated, and m and n are two ages. Thus, the null hypotheses are β = 0 and λ = 0, which reflect a stationary covariance function.

In our analyses, we used the expected covariance matrix for a given pedigree [3], which has elements giving the covariance in trait value *y* for any two relatives:

where x and z index individuals at ages m and n, respectively, Φ_{
xz
}gives their kinship status, π is the probability that a random allele is IBD at the QTL, δ_{xz} is 1 for x = z and 0 otherwise, subscripts g, q, and e denote the polygenic, QTL, and environmental components, respectively, and average age is taken over the sample population, defined as all individuals measured for the trait of interest. Some points of clarification are needed. This is a cross-sectional model that applies generally to three types of pair-wise comparisons of individuals. In one type, let x = z such that *m* = *n*. Equation (5) gives the variances in this situation, in accord with the standard variance components model. In a second type, it may be such that x ≠ z while m ≠ n, and, in a third type, it may be such that x ≠ z while m ≠ n. Note that none of these types are longitudinal comparisons, which would be the case where x = z while m ≠ n (i.e., the same individual is measured at different ages). The goal of this approach lies in estimation of the change parameters, namely β and λ, so that we can test the null hypotheses stated above.

For the bivariate model, a trait measured at two different time points is treated as a bivariate phenotype. The pedigree covariance matrix model for a bivariate phenotype K, with constituent phenotypes I and J, may be written as [4]:

**Σ**_{
K
} **= 2Φ****G +** **Q + I****E,** (6)

where the new matrices **G**, **Q**, and **E** convey the polygenic, QTL, and environmental variance components, respectively, and is the Kronecker-product operator. The variance components for this model are those for the univariate model for traits I and J, given by Σ_{I} and Σ_{J} (cf. equation (1)), and the trait cross-covariances, which may be parameterized as ρ_{gij}σ_{gi}σ_{gj}, ρ_{qij}σ_{qi}σ_{qj}, and ρ_{eij}σ_{ei}σ_{ej} for the polygenic, QTL, and environmental components, respectively [4]. The null hypotheses of no polygenic genotype or QTL × age interaction are expressed as ρ_{gij} = 1 and ρ_{qij} = 1, respectively.

Scatter plots were examined for traits showing increasing or decreasing variance in trait values with age, suggesting potential for G × age interaction. Traits meeting these criteria, namely systolic blood pressure (SBP) and fasting glucose (GLUC), were analyzed using SOLAR [2], with age, sex, hypertension medication, and body mass index (BMI) as covariates. Phase 1 analyses were conducted on an augmented data set. Exams 12 (1970), 16 (1978), 18 (1982), and 20 (1986) in Cohort 1 were combined with Exams 1 (1971), 2 (1979), 3 (1983), and 4 (1987) in Cohort 2, respectively. To reduce kurtosis after combining, a few outliers were removed for SBP (1–6 per exam, all with values > the mean) but many for GLUC (51–101 per exam, all but one > the mean). No data transformation was made. Since diabetic status was not available, we could not control for it. This gave combined measurement periods 1–4.

Based on genome scans that we conducted for both SBP and GLUC, we decided to focus on SBP. For Phase 2, data for SBP values corrected for hypertension treatment for Cohorts 1 and 2, imputed following Levy et al. [9] (see Soler and Blangero [10]), were analyzed. Bivariate analyses were performed on residuals from multiple regressions with corrected SBP values as the dependent variable and age, sex and BMI as independent variables. For the bivariate analyses, analysis of residuals was necessary due to the unbalanced and incomplete structure of the FHS longitudinal data, which made covariate specification unrepresentative across measurements. For this SBP data set (corrected values and residuals), we experimented with combined exams from Cohorts 1 and 2 and with exams from Cohort 2 taken separately.