Data
We analyzed the longitudinal Framingham Heart Study Family Cohort data set, provided for the Genetic Analysis Workshop, which included 4692 individuals in all 330 pedigrees. The primary response variable of the GLMMs and linkage models reported in this paper was the SBP measured over time. Explanatory covariates for all models included: sex, age, height, weight, number of drinks per week, number of cigarettes per day. Sex (m = 1, f = 0) was analyzed as a binary covariate. Other variables were treated as continuous. All continuous covariates were centered at or close to their mean. In order to allow for a nonlinear increase in SBP with age, quadratic and cubic age terms were also included. Although the model is Bayesian, the effects of each covariate on the phenotype can also be assessed using a pseudo-likelihood approach.
The covariates drink (number of drinks per week) and CPD (number of cigarettes per day) were only measured at 11 and 18 time points, respectively, in the first cohort. However, the majority of individuals had only one or two distinct values for these covariates throughout the study, so each missing value of drink and CPD was replaced with the mean of the covariates for that individual in the two observations surrounding the missing observation. This enabled us to increase the amount of useful observations without excluding important covariates from the model. The phenotypic component of the model included 24,746 observations for 2860 individuals. These observations all had non-missing values for SBP, height, weight, age, drink, and CPD.
Variance components model
The models fitted in this analysis are extensions of our previously described variance components models for extended pedigrees with normal, binary, or censored survival phenotypes [4, 5]. The original models are GLMMs with components of variance due to additive polygenic effects (σ2A), common family environment (σ2C), and common sibling environment (σ2Cs), and since the phenotype in this analysis is assumed normally distributed, an error variance (σ2E). The complex familial correlation structure is modeled through the use of shared random effects for each variance component. The model for the jth individual in the ith pedigree can be represented as:
μij = β0 + βTXij + Aij + Cij + Csij
yij ~ N(μij, σ2E),
where yij is the observed phenotype, Xij is a matrix of covariates, Aij ~ N(0, VA), Cij ~ N(0, VC), Csij ~ N(0, VCs) and the random effects are sampled from appropriate multivariate normal distributions with variance-covariance matrices (VA, VC, and VCs) parameterized by σ2A, σ2C, and σ2Cs, respectively, and structured as described in [5]. Individuals without phenotypic data are used in generating the random effects but not in the main model.
The extensions not only allow for repeated measures on individuals, but also use the extra information available to assess evidence for genetic effects on the rate of change of the outcome over time [6]. At the kth observation for the jth individual in the ith pedigree, the model can be expressed as:
μijk = β0 + βTXijk + Aij + Cij + Csij + ageijk × (βage + A.timeij + C.timeij + Cs.timeij + E.timeij)
yijk ~ N(μijk, σ2E),
where yijk is the observed phenotype, Xijk is a matrix of covariates, βage is the mean gradient of age (for all individuals), Aij ~ N(0, σ2A), Cij ~ N(0, σ2C), Csij ~ N(0, σ2Cs), A.timeij ~ N(0, σ2A.time), C.timeij ~ N(0, σ2C.time), Cs.timeij ~ N(0, σ2Cs.time), E.timeij ~ N(0, σ2E.time).
In this model time was measured using the covariate age. The random effects for the gradient of age, A.timeij, C.timeij, Cs.timeij, and E.timeij, have a covariance structure equivalent to that described previously [5], with each variance component replaced by its counterpart for the gradient, so (for example) in this part of the model σ2A is represented by σ2A.time. However, the variance components for the intercepts (σ2A, σ2C, σ2Cs, and σ2E) are independent of the variance components for the corresponding gradients (σ2A.time, σ2C.time, σ2Cs.time, and σ2E.time), and the estimated random effects for each variance component for the intercept are different from the estimated random effects for each variance component for the gradient. All random effects are constant across observations on the same individual; they do not change over time. Just as the size of σ2A can indicate whether additive genetic effects influence the outcome, the size of σ2A.time can indicate whether additive genetic effects influence the rate of change of the outcome with time. The narrow-sense heritability (h2N) was defined as the ratio of variance due to additive genetic effects (σ2A or σ2A.time) to the total phenotypic variance in SBP or rate of change in SBP [7].
The models were fitted using Gibbs sampling in WinBUGS v1.3 [8]. Because WinBUGS uses a Bayesian formulation, prior distributions must be specified for all parameters of interest. Vague normal N(0, 1000) priors were specified for fixed regression coefficients. Vague Pareto(1, 0.001) priors were specified on the precisions (inverses of the variances) of all random effects; these were equivalent to specifying uniform priors on the corresponding variances. Models were run for 45,000 iterations after a burn-in of 5000 iterations.
Blood pressure adjustment
The problem of how to model continuous SBP when some individuals are on blood pressure treatment is complex. However, since complex modelling of SBP was not the main aim of this particular analysis, we chose a simple method of adjustment based on known average treatment effects [9, 10]. This involved adding a constant (10 mm Hg) to each phenotype where the individual was on treatment, to reflect the 'true' SBP that might have been observed if the individual had not been on treatment.
Linkage analysis
As described previously, the random effects due to σ2A (Aij, the sigma-squared A residuals, or SSARs) can be used as a continuous phenotype in a linkage analysis [11, 12]. The SSARs for the gradient (A.timeij, the sigma-squared A time residuals or SSATRs), may also be used in such a way. The estimated SSARs and SSATRs were used as continuous phenotypes in a multipoint variance components linkage analysis that was performed using SOLAR v1.7.4 [13]. All 395 markers were included. The most common approach undertaken in previous genetic analyses of longitudinal data has been to estimate a simple linear gradient in each subject and then use this as a phenotype. Therefore, in order to compare our SSAR methodology to such an approach, we also estimated a simple linear gradient in each subject and used this as an outcome into the SOLAR linkage analyses.