Abstract
Longitudinal studies are an important tool for analysing traits that change over time, depending on individual characteristics and environmental exposures. Complex quantitative traits, such as lung function, may change over time and appear to depend on genetic and environmental factors, as well as on potential geneenvironment interactions. There is a growing interest in modelling both marginal genetic effects and geneenvironment interactions. In an admixed population, the use of traditional statistical models may fail to adjust for confounding by ethnicity, leading to bias in the genetic effect estimates. A variety of methods have been developed to account for the genetic substructure of human populations. Familybased designs provide an important resource for avoiding confounding due to admixture. To date, however, most genetic analyses have been applied to crosssectional designs. In this paper, we propose a methodology which aims to improve the assessment of main genetic effect and geneenvironment interaction effects by combining the advantages of both longitudinal studies for continuous phenotypes, and the familybased designs. This approach is based on an extension of ordinary linear mixed models for quantitative phenotypes, which incorporates information from a caseparent design. Our results indicate that use of this method allows both main genetic and geneenvironment interaction effects to be estimated without bias, even in the presence of population substructure.
Keywords:
geneenvironment interaction; longitudinal phenotypes; power; bias; population substructureIntroduction
In spite of the multiple efforts to find genetic factors conferring susceptibility to complex diseases, the success of genetic association studies is still hampered by the difficulty in replicating findings in different populations. Among the plausible explanations for this lack of replication is the fact that the effects of environmental factors, which can interact with genetic factors, are not always taken into consideration [1]. There is an increasing interest in studying different susceptibilities to environmental factors in subjects with different genotypes; however, power and bias issues with regard to the statistical estimation of geneenvironment interaction effects persist.
Highquality information about individual environmental exposure is crucial for the assessment of geneenvironment interactions [2]. Failure to measure changes in exposure levels over time could lead to an underestimation of the role of the environment in the interaction. Repeated measurements of the temporal relationship between an outcome and the exposure may overcome such a problem when both the endpoint and the exposure are timedependent variables. In addition, potential misclassification due to ambiguity in the definition of complex diseases may be avoided through the measurement of quantitative diseaserelated phenotypes as the outcomes of interest. For example, quantifying the decrements in lung function over time through repeated spirometric tests may provide insights into the pathogenesis of chronic obstructive pulmonary disease (COPD) or asthma. Many disease 'predictor' phenotypes are thought to change withinsubject because of both environmental and genetic factors, and of their potential interactions over time.
On the genetic side, population substructure is an important practical issue for genetic association studies. When the study population is not a collection of randomly mating individuals, several discrete subgroups that are genetically different may be identified; the collection of these subpopulations is referred to as population substructure or stratification [3]. Moreover, disease prevalence also tends to differ among these subgroups [4]. Consequently, without stratification adjustment, allele frequency can appear to be associated with the disease, regardless of whether the genotype has a functional effect on that health outcome or not. By contrast, when the genotype distribution is homogeneous among groups, population substructure may not be an issue. For example, if people are randomly assigned to treatment groups, it is expected that those groups will be genetically similar. If, additionally, there are no differences in the response to treatment among the different subgroups, bias due to population substructure is unlikely.
Another source of spurious associations is population admixture, which refers to the mixture of different ancestries; that is, people from different ethnic groups interbreed, so the genome of the new generations is a combination of genotypes of the original ancestry groups, and, consequently, in some genes, allele frequencies are not homogeneously distributed in the study population. For example, it has been recognised that Latino populations have varying proportions of African, Native American and European ancestry [5]. Like population substructure, if the risk of disease depends on ancestry, a high risk of disease may be erroneously associated with a high allele frequency; thus, in admixed populations, ethnicity may confound associations between genotype and outcome and assessment of gene by environment interactions. The direction of the confounding could be positive or negative. Therefore, to identify true associations, population substructure must be taken into account in the analysis.
With the increasing availability of genetic data, there is a growing interest in modelling both marginal genetic effects and geneenvironment interactions. Inclusion of interactions, when they exist, can increase the statistical power of detecting both genetic and environmental effects [6]. Traditional statistical models for detecting significant main effects and interactions may not be completely adequate for studying genetics in admixed or stratified populations, however.
A variety of methods have been developed to account for the genetic substructure of human populations [7]. Familybased designs provide an important resource for avoiding confounding due to admixture [8]. The simplest design for testing association is the caseparent (or trio) design because it uses genotypes from an affected offspring, the case, and his/her two parents. The outcome is measured, however, only in the offspring. Many of these methods have been developed for crosssectional designs, but can be applied to repeated measurements through the twostep modelling approach. The first step consists of calculating the slope between the longitudinal outcome and the timedependent environmental exposure; thus, we calculate a single individual endpoint, the slope, for each subject. In the second step, the genetic methods for crosssectional studies, where the slope is the single outcome, can be applied [9].
In this paper, we first provide a short review of different approaches for studying geneenvironment interactions for quantitative traits, and then propose a method that aims to improve the assessment of main and geneenvironment interaction effects by combining the advantages of both longitudinal studies for continuous phenotypes and the familybased designs. This approach is based on an extension of ordinary linear mixed models (OLMM) for quantitative phenotypes which incorporates information from a caseparent design. We call the model the 'adjusted linear mixed model' (ALMM), and through simulation methods we show that even when population stratification is present, both main genetic and geneenvironment interaction effects can be estimated without bias, and that this is more powerful than the twostep modelling approach.
The broad objectives of this paper do not extend to giving technical details about the familybased approach and its extensions, or to giving an extensive explanation about linear mixed models. Rather, we present what we consider to be a widely applicable method for correctly assessing the main genetic effect and geneenvironment interactions for timedependent quantitative traits in stratified populations. For this purpose, we use simulated repeated measurements of forced expiratory flow between 75 per cent and 25 per cent of vital capacity (FEF_{2575}) ie (lung function) on asthmatic children exposed to ozone pollution, based on the observed distributions in a real cohort study conducted in Mexico City [10].
In order to set the stage for our methodology, we first provide a brief overview of some existing ordinary linear regression (OLR) models for testing main genetic effects and geneenvironment interactions in crosssectional studies that incorporate information about parental genotype (caseparent or trio design), adjusting for admixture. We then briefly present the familybased association test (FBAT) approach, which, as a second step (after computing the slope between the outcome and the exposure), represents an alternative method for analysing genetic associations over time. We next review the ordinary linear mixed models (OLMM) which are a standard approach for the analysis of longitudinal data, and present the adjusted linear mixed models (ALMM) as an extension of OLMM combined with the adjusted crosssectional regression models. In order to show that the twostep modelling approach provides a valuable alternative for analysing longitudinal data, we explain the relationship between this approach and the linear mixed models. Finally, we give details about the simulation procedures and present our results and discussion.
Methods
Models for crosssectional data with a single quantitative measure for each subject
Existing methods for testing main genetic effect and geneenvironment interactions with a single measured outcome include (1) OLR models and extensions that aim to adjust for ethnicity by including parental genotype information in a caseparent or trio design, and (2) the FBAT approach, which uses a score test based on a conditional likelihood [11].
The OLR approach
In a genetic association analysis with quantitative traits that follow a linear model, the assessment of genetreatment interactions may be conducted using standard linear regression models for independent subjects. Under the usual assumptions, the well known model for testing the interaction between two covariates is:
where:
n is the number of subjects in the study;
X_{i }is a fixed variable that translates an offspring genotype into a numeric value; and
Z_{i }is an observed environmental covariate, either continuous or dichotomous.
X_{i }is a scalar whose value depends on the disease model. If the locus has two alleles, A and a, the additive model assumes that each copy of the variant allele* 'A' changes the outcome in an additive amount. Thus, X_{i }counts the number of 'A' alleles in the offspring genotype (X_{i }= {0, 1, 2}). In the recessive model, X_{i }is coded as an indicator variable for the AA genotype. As a special case, model (2) is used for testing main genetic effects adjusted for the environmental covariate:
A rejection of the null hypothesis H_{0 }: β_{1 }= 0 means that the quantitative trait is associated with the alleles in the marker.
The caseparent or trio design
Unlike OLR models, family designs aim to avoid spurious associations due to population admixture. In the caseparent design, the proband is the offspring that identifies the family for the study; the genotypes of the candidate gene are measured for all members of the trio, but the quantitative trait is measured only in the offspring. The alternative form of the OLR approach for testing main genetic effect on quantitative outcomes, where the parental genotype information is included, was developed by Allison;[12] it is based on the simplest familybased design for testing associations, known as the transmission disequilibrium test (TDT) [13]. The model is adjusted for the expected value of the offspring's genotype conditional on the parental genotypes; thus, the adjusted version of model (2) is:
where:
g_{im}, g_{if }are the parental genotypes (mother and father, respectively) and E(X_{i}g_{im}, g_{if}) is calculated under segregation and independent assortment assumptions using Mendel's law. Its value depends on the mating type and the disease model.
The adjusted genotype X_{i } E(X_{i}g_{im}, g_{if}) is the subject's deviation from the family mean under Mendel's law. β_{1 }represents the withinfamily effect of the gene on the outcome. As a result of centring X_{i }by its expected value conditional on parental genotypes (g_{im}, g_{if}), ethnicity bias is avoided, since all possible genotypes  depending on the mating type  are taken into account, even those that were not transmitted to the affected offspring. This strategy does not, however, necessarily prevent bias due to other kinds of population stratification,[14] such as the one that occurs when parental mating type is highly correlated with the levels of exposure, for example. For this reason, Allison[12] and Ewens et al.[14] propose an alternative version of model (3) in which the intercept, representing the betweenfamily component, depends on the mating type as a fixed effect:
where:
and M = 1, 2, 3 are the three possible mating types with at least one heterozygous parent, including informative families only.^{† }Note that here and in subsequent equations, M depends upon i via the parental genotypes, but this is suppressed for simplicity.
It is noteworthy that since both β_{0M }and E(Xg_{im}, g_{if}) are constant within the mating type, the estimation of the main genetic effect (β_{1}) through model (4) is completely equivalent to using model (5):
Following the same idea, Gauderman[15] proposed a likelihood ratio test (LRT) of β_{1 }in model (5), although, in order to increase the power for detecting main genetic effects, he included the whole set of families, regardless of the heterozygous condition. This is called the quantitative transmission disequilibrium test with mating type indicators (QTDT_{M}). An extension of the QTDT_{M }to include the environmental covariate and the geneenvironment interaction was also suggested by Gauderman[15] with the model:
Where M = 1, 2,..., 6 are the six possible matingtypes.
In contrast to the relationship between models (4) and (5), model (6) is not equivalent to:
because, although β_{0M }and E(X_{i}g_{im}, g_{if}) are constant within mating type, the environmental covariate (Z_{i}) is not.
An alternative model to (7), (and an extension of Gauderman's idea), is given by the adjusted model that we call the adjusted quantitative transmission disequilibrium test with matingtype indicators (AQTDT_{M}).
where now both the intercept and the slope for the environmental covariate depend on the mating type, which means that the model has been adjusted for a potential correlation between the exposure to environmental factors and the mating type. Note that now model (8) is equivalent to:
The advantage of using models (5) and (9) is that the inclusion of an indicator variable for the mating type, rather than calculating the expected genotype, provides protection against population admixture while taking into account those situations where environmental exposure (Z) may depend on the mating type and other additional types of population substructure [14]. As usual, in (6), (7), (8) and (9), β_{3 }estimates the effect modification of the gene on the environmental effect Z_{i }and H_{0 }: β_{3 }= 0 is the null hypothesis that states the no interaction effect.
The FBAT approach
FBAT is the generalisation of the TDT for the trio design. It encompasses a broad class of statistical methods for testing genetic associations adjusting for potential admixture or stratification. Such methods are also based on extensions of the TDT and regression models, although the covariance between genotype and phenotype is the statistic of interest [11]. The general FBAT statistic has been explained elsewhere [16]. Briefly, for n nuclear families, one offspring in the family i and no covariates:
where:
and E(X_{i}g_{im}, g_{if}) and Var(X_{i}g_{im}, g_{if}) are calculated under the null hypothesis of Mendel's law. This statistic follows a chisquare distribution with one degree of freedom. In addition, unlike regression models where the offspring's genotype is assumed to be fixed and observed, the general FBAT approach considers the offspring genotype as a random variable. The general idea is first to calculate a test statistic for the association between the trait and the marker locus, and then, as a second step, the distribution of this test statistic is derived from the distribution of the offspring genotype under the null hypothesis of no association. The distribution of the test statistic is computed conditioning on the sufficient statistic given by the parental genotype and the observed offspring's phenotype. Under these conditions, no assumptions about the allele frequency, the recombination fraction or the penetrance function are required. Due to the fact that the general FBAT statistic can only test main genetic effects, and since the test statistic is based the relative size of U with respect to its standard deviation, the genetic effect is not directly estimated.
Under the philosophy of the FBAT, Vansteelandt et al.[17] proposed an extension that permits the assessment of and testing for the geneenvironment interaction without any assumptions about the genotype distribution and with no bias due to unmeasured ethnicity confounding through which Mendel's laws hold. Such an extension is based on Gestimation (causal inference) methodology and is called QBATI; this test statistic has the same form as the general FBAT (10), although the expected genotype and the U statistic have more complex expressions. A brief presentation of both FBAT and QBATI can be found in Appendix 1.
FBAT and QBATI statistics are available from PBAT free software at http://biosun1.harvard.edu/~clange/pbat.htm webcite and in the library pbatR under the R package environment.
Models for longitudinal data
Until now, we have discussed methods for analysing geneenvironment interactions when a single measurement of the phenotype is available. We now turn to a discussion of methods for the analysis of quantitative repeated phenotype data to evaluate the effects of a gene and the environment over time.
In longitudinal designs, the unit of study is not each individual or each measurement, but rather the sequence of measurements on each subject. This means that the major advantage of a longitudinal analysis is that the socalled cohort and age effects are estimated separately; that is, differences among people in their baseline levels (cohort effects), can be discriminated from the changes over time (ageing effects) within individuals. In other words, measurements across people and repeated values across time are sources of strength. Note that crosssectional data provide information for assessing only the former effect; thus, longitudinal studies tend to be more powerful than crosssectional studies [18].
Although there are different approaches to longitudinal data analysis, we consider here the two most commonly used: OLMM and twostep methods.
OLMM
OLMM assumes that the vector of repeated measurements on each subject follows a linear regression model. Thus, each individual model may have subjectspecific intercept and slope (the random effects) representing the different susceptibilities to the environmental exposure among subjects. For most outcomes, variability across individuals is greater than withinsubject. This difference may be due the influence of genetic composition [19].
In general, a linear mixedeffects model satisfies:
where:
i = 1; 2; 3;..., n subjects;
j = 1; 2; 3;..., m corresponding times at which the measurements are taken on each subject;
t_{ij }is the repeated time (or exposure) variable;
b_{1i }is the random subject intercept effect (α_{0 }+ b_{1i}), which varies among subjects;
b_{2i }is the random subject slope effect (α_{1 }+ b_{2i})t_{ij}, which varies among subjects;
e_{ij }is a random variable regarded as measurement or sampling errors;
and X_{i }and Z_{i }are as previously defined. Note that model (11) assumes that both genotype (X_{i}) and environmental exposure (Z_{i}) remain fixed over time. However, we can take advantage of the time variable (t_{ij}) to include an extra source of exposure that changes across time. For example, we are particularly interested in the assessment of the geneenvironment interaction effect between the glutathioneStransferase M1 gene (GSTM1) and antioxidant supplementation on lung function  FEF_{2575 } of asthmatic children exposed to ozone.
Therefore, in model (11), X_{i }represents the individual GSTM1 genotype and Z_{i }is a dichotomous variable that denotes the antioxidant supplementation group which was randomly assigned at baseline and remained fixed during followup, and which since ozone is a timedependent variable, can be represented by t_{ij}. This is a model with random intercept and random slope on ozone. As for ordinary linear regression models in crosssectional studies, OLMM assumes independent subjects and does not account for population admixture. Therefore, in order to avoid potential estimation bias, it is necessary to adjust for ethnicity.
ALMM
The ALMM form is a straightforward extension of the approaches presented for the crosssectional designs. That is, following the ideas of Allison[12], Ewens et al.[14] and Gauderman[15], model (12) can be rewritten using the indicator variables for the mating type and the offspring's conditional expected genotype:
where i, j, b_{1i}, b_{2i}, e_{ij }are as defined for (11); X_{i }represents the GSTM1 genotype; Z_{i }denotes the antioxidant supplementation group; t_{ij }represents the ozone exposure and M = {1, 2, ..., 6} are the six possible mating types. This model is an extension of (11) and is able to control for population admixture and for any dependence of the environmental exposure on mating type. Once again, using an indicator variable for the mating type prevents controls for other potential sources of population structure.
Note that, as in model (8), the simplest, but equivalent, expression for the above model that does not need the calculation of the expected genotype is given by:
Through this model, the main genetic and environmental effects and geneenvironment interactions can be assessed using standard statistical software. Additional advantages of this model are that genetic and longitudinal effects are estimated simultaneously; thus, parameters are mutually adjusted for each other[9] and within and betweenindividual variabilities are taken into account. Moreover, the model can be adjusted by other covariates, such as gender and age at baseline, or even by smooth functions of variables which have nonlinear association with the outcome. In particular, in the example, two sources of environmental exposures  supplementation treatment (a constant) and ozone (which is time dependent)  are included. Computationally, longitudinal models require more complex algorithms than those corresponding to crosssectional designs; however, this is no longer a problem for the users of statistical packages. Appendix 1 includes the R code for estimating the effects of model (13).
Twostep modelling approach
As the term implies, this approach includes two separate steps. In the first step, it is assumed that the repeated measures on subject i, are independent and follow a linear regression model. Therefore, the individual intercept and slope are estimated by an ordinary linear regression model:
Therefore, different regression coefficients (γ_{0 }intercept and γ_{1 }slope) correspond to different individuals. For example, with the lung function data, the outcome is the repeated FEF_{2575 }and γ_{1 }represents individual susceptibility to ozone over time. Thus, the longitudinal observations are reduced to one summary statistic per subject [20]. The second step includes the genetic analysis where these slopes are the single outcome for each person. This approach has been frequently used in segregation, linkage and association analysis [9]. The R code for these models is included in Appendix 1.
It is interesting to note the close relationship between the coefficients in the twostep model and the respective coefficients in a linear mixed model. That is, by definition, from the ordinary linear regression model (14):
where Y_{ij }is given by (11) and
Thus, it is straightforward to show that
This therefore enables one to relate the α coefficients of the linear mixed model to the β coefficients of the ordinary linear regression model given by (1) as follows: β_{1 }= α_{5 }and β_{3 }= α_{7 }representing the main genetic effect and the genetreatment interaction on the slope, respectively.
This approach is remarkably similar to that of mixed models, although longitudinal and genetic effects are not jointly estimated, and timedependent covariates need to be summarised in one measurement  the mean or the median, for example  in order to be included in the analysis. Computationally, this procedure is simpler than mixed models, and any elemental statistical package will suffice for conducting the analysis.
Table 1 summarises the different models included in this section.
Table 1. Regression models included in the paper. The 'Model' column refers to the number that identifies each model in the paper.
Simulations
Here, we address differences in the analytical approaches, in terms of both bias and power, for detecting main genetic effects and geneenvironment interactions using simulations.
In order to compare the different methods presented, we simulated data with similar characteristics to those in the study by Romieu et al.[10] Briefly, this study was a randomised trial using a doubleblinded and longitudinal design, including antioxidant supplementation for asthmatic children who were residents of Mexico City and therefore exposed to ozone pollution. There were 12 repeated measures for both FEF_{2575 }and ozone. The deletion polymorphism of GSTM1, absent versus present, was determined for each child and, through a stratified analysis, evidence for interaction between the antioxidant treatment (dichotomous and fixed variable) and the GSTM1 genotype was seen for the effect of ozone on lung function.
For our simulations, 12 repeated measures of lung function (FEF_{25  75}) were generated for each offspring with a mean vector given by a linear mixed model conditional on treatment, genotype and ozone level, assuming an additive (or recessive) disease model (15). An error vector was added to each mean by drawing from a multivariate normal distribution. The variancecovariance matrix of the errors was assumed to be equal to the observed variancecovariance matrix among the residuals in the real data, where model (15) was fit to the repeated FEF_{25  75 }measurements.
For the purpose of using the familybased approach, samples of independent trios were simulated. Each parent was randomly assigned a genotype assuming the HardyWeinberg equilibrium, while each offspring was assigned a genotype assuming Mendel's law. Treatment (Z) was randomly and independently assigned for each subject i with a 50 per cent probability for supplement or placebo group. Both additive and recessive disease models were considered.
With regard to the population stratification, two different situations were considered: the first one assumed a homogeneous population (HP), where the observed allele frequency for GSTM1, P(a) = 0:4 and P(A) = 0.6,^{‡ }was used for simulating the genotypes. The generating model we used had the following form:
Where: t_{ij }represents the timedependent ozone exposure; Z_{i }is a dichotomous variable representing the treatment group and X_{i }is either a continuous variable counting the copy variant number of GSTM1 (0, 1 or 2) in the additive disease model or a dummy variable in the recessive case. The parameters are:
α_{0}, α_{1}, α_{2}, α_{3}, α_{4 }are equal to the corresponding estimates obtained through a linear mixed model with random intercept and random slope on ozone applied to the original dataset. With the only aim to simulate situations where the effects are clearly identified without using large samples but rather looking at the differences between the methods, α_{5}, α_{6}, α_{7 }were magnified. Assigned values to α_{5}, α_{6}, (α_{7}) try to mimic the respective effects found in the stratified analysis reported by Romieu et al.[10] Specifically, the genotype effect for each part per billion (ppb) of ozone was 0.8 ml/s in the placebo group (n = 78 subjects) and 0.16 ml/s in the supplement group (n = 80); thus, the main genetic effect was taken as the rounded weighted average effect, while ignoring the supplement effect. The same procedure was applied to calculating the main treatment effect. Finally, since the effect of antioxidants was stronger in the GSTM1 null genotype group (0.95 ml/s), and there was no significant effect in the GSTM1positive group, the coefficient for the genetreatment interaction on the lung functionozone relationship was rounded to 1.0. Table S1 (Table 6) summarises the models used in the simulation process and Table S2 (Table 7) shows the observed effects in the real cohort study conducted in Mexico City.
Table S1. List of models used in the simulation process. The column model refers to the number that identifies each model in the paper. X_{i }is a fixed variable that translates an offspring genotype to a numeric value; Z_{i }is an observed environmental covariate, either continuous or dichotomous; g_{im}, g_{if }are the parental genotypes (mother and father, respectively); E(X_{i}g_{im}, g_{if}) is calculated under segregation and independent assortment assumptions using Mendel's law; M = 1, 2, ..., 6 are the six possible mating types; i = 1, 2, 3, ..., n subjects; j = 1, 2, 3, ..., m measurement occasions into the subject; t_{ij }is the repeated ozone exposure variable.
Table S2. Observed effects in the real cohort study conducted in Mexico City. 95% CI = 95% confidence interval.
The second situation presumes a stratified population (AP) with a 50/50 mix coming from two populations with different allele frequencies and different susceptibilities: P_{1}(A) = 0.4; P_{1}(a) = 0.6 and P_{2}(A) = 0.8; P_{2}(a) = 0.2. Note that, on average, the combined population has the same allele frequencies as the homogeneous one. With the purpose of simulating differential susceptibilities to ozone and supplementation, although allowing the bias assessment for the main genetic and interaction effects, generating model (15) differs in β_{0 }and β_{6 }coefficients (based on the observed percentiles). That is, in the first population, the observed 95th percentile for FEF_{25  75 }(α_{0 }= 3.3) was used as the intercept, and no treatment effect on the slope α_{6 }= 0 was assumed. On the other hand, in the second population, the 5th percentile (α_{0 }= 0.75) was taken as the intercept, and a strong treatment effect on the slope was assumed, α_{6 }= 2 (meaning that 20 ml/s/10 ppb decreased lung function, on average, in the placebo group in comparison with the supplement group). The variancecovariance matrix was constant over the different simulated samples.
Regarding the genetic effect, two different scenarios were considered. The first scenario represents situations where the variability in the outcome may be attributable just to the main effect of the gene and the treatment, meaning that there is no genetreatment interaction effect (α_{7 }= 0) in the generating model (15). The second scenario assumes that all genetic, treatment and genetreatment interaction are present in the true model.
Assessment of main genetic effect
The first scenario, where there is no genetreatment interaction (α_{7 }= 0), was used for testing the main genetic effect, adjusted by treatment effect.
Under the twostep modelling approach, the slope between the outcome and ozone was first computed. In the second step, the ordinary linear regression model (16), the AQTDT_{M }(17) model and the FBAT statistic (18) were used for testing the null hypothesis of no main genetic effect (H_{0}: β_{1 }= 0).
The corresponding null hypothesis when using longitudinal outcomes H_{0}: α_{5 }= 0 was tested through ordinary and adjusted linear mixed models; that is:
and
For the purpose of comparison, statistics based on model (16) will be referred to as OLR, those based on (17) will be referred to as the QTDT_{M}, those based on (19) will be referred as OLMM and those based on (20) will be referred as ALMM.
Assessment of geneenvironment interaction
Using the same idea, the second scenario, where α_{7 }= 1 in the generating model, was used for assessing the interaction effect through a one degree of freedom test.
Assuming one outcome per individual, in the twostep modelling approach, the null hypothesis H_{0}: β_{3 }= 0, was tested using (21), (22) and the QBATI statistic:
For repeated measurements, the corresponding null hypothesis H_{0}: α_{7 }= 0 was tested using the OLMM and ALMM models respectively:
and
For the purpose of comparison, statistics based on model (21) will be referred to as OLR, those based on (22) will be referred to as the adjusted AQTDTM, those based on (23) will be referred as OLMM and those based on (24) will be referred as ALMM.
The empirical power for each test was estimated as the percentage of occasions on which the null hypothesis was rejected at a significance level α ≤ 0.05 for a twosided test. In each simulation study, 1,000 independent replicate datasets were generated. Each dataset consisted of n (n = 100, 200, 300, 400, 500, 600 and 1,000) complete and independent trios.
Bias was calculated as the average of the difference between the estimator and the true parameter value .
Results
Estimation bias
In order to look at the differences among methods, the estimation bias for the main genetic effects and geneenvironment estimation was computed under both population conditions, homogeneous (HP) and stratified (SP) populations. Table 2 shows the resultant bias for the four methods, with two columns per method (HP and SP). When there is no ethnicity confounding (HP), all methods (OLR, AQTDT_{M}, OLMM and ALMM) for estimating main effects are unbiased, regardless of the design analysis (independent subjects or trios) or the modelling approach (twostep or longitudinal data). By contrast, when the population is stratified, the selection of the design is crucial. In other words, while estimators obtained from the caseparent design are unbiased, models using independent subjects underestimate the effects by around 0.36 units. Note that, regardless of the modelling approach, both designs provide similar results.
Table 2. Bias results for main genetic effect assessment comparing ordinary statistical methods (OLR and OLMM) to familybased methods (AQTDT_{M }and ALMM) under homogeneous (HP) and stratified (SP) populations. Each time, n cases were simulated with parameters β_{1 }= α_{5 }= 0.5. Simulations are based on the additive genetic model. † = number that identifies each model in the paper.
As in the main genetic effect assessment, the estimation of the geneenvironment interaction also depends strongly on the design of the study when there is population stratification. In this case, when using ordinary statistical methods, the interaction effect was also underestimated by about 0.72 units, but with a homogeneous population the design is not relevant in terms of bias (Table 3).
Table 3. Bias results for geneenvironment interaction effect assessment comparing ordinary statistical methods (OLR and OLMM) to familybased methods (AQTDT_{M }and ALMM) under homogeneous (HP) and stratified (SP) populations. Each time, n cases were simulated with parameters β_{3 }= α_{7 }= 1. Simulations are based on the additive genetic model. † = number that identifies each model in the paper.
Empirical power
Since power comparisons among biased and unbiased methods cannot be fair, the power of OLR and OLMM under population stratification was not computed. Regarding genetically homogeneous populations, in both main genetic and genetreatment interaction effects, ordinary regression models are the most powerful methods (Figures 1 and 2). Moreover, and as was expected, the use of repeated measures (OLMM) is more powerful than the use of a single measure (OLR). Note that, among the familybased models, ALMM is the most powerful and that FBAT or QBATI statistics and AQTDT_{M }are equivalent with regard to power.
Figure 1. (a). Empirical power results for main genetic effect assessment comparing different methods under the assumption of a homogeneous population. Simulations are based on the sample size as indicated in the plot and the additive genetic model. (b) Empirical power results for geneenvironment interaction effect comparing different methods under the assumption of a homogeneous population. Simulations are based on the sample size as indicated in the plot and the additive genetic model.
Figure 2. (a). Empirical power results for the main genetic effect assessment comparing different methods under the assumption of a stratified population. Simulations are based on the sample size as indicated in the plot and the additive genetic model. (b). Empirical power results for the geneenvironment interaction effect comparing different methods under the assumption of a stratified population. Simulations are based on the sample size as indicated in the plot and the additive genetic model.
When the population is genetically mixed, all methods, regardless of the modelling approach or the design, lose power in comparison with the setting where the population is homogeneous (Tables 4 and 5). Once again, ALMM is the most powerful method for detecting main genetic or the interaction effects. In other words, the modelling approach may be crucial, in terms of power, for conducting the association analysis.
Table 4. Empirical power results for main genetic effect assessment comparing ordinary statistical methods (OLR and OLMM) to familybased methods (AQTDT_{M}, FBAT and ALMM) under homogeneous (HP) and stratified (SP) populations. Each time, n cases were simulated with parameters β_{1 }= α_{5 }= 0.5. Simulations are based on the additive genetic model. † = number that identifies each model in the paper.
Table 5. Empirical power results for geneenvironment interaction effect assessment comparing ordinary statistical methods (OLR and OLMM) to familybased methods (AQTDT_{M}, QBATI and ALMM) under homogeneous (HP) and stratified populations (SP). Each time, n cases were simulated with parameters β_{3 }= α_{7 }= 1. Simulations are based on the additive genetic model. † = number that identifies each model in the paper.
Additive versus recessive disease models
In the recessive model, the same relationship among methods is observed, although the additive model is always more powerful than the recessive one, regardless of the testing approach. The data are shown in the supplementary tables. This may be related to the fact that in the additive model, X has a wider range of variation, whereas in the recessive one, X is an indicator variable. In addition, the number of informative families for FBAT statistics is always larger when the additive disease model is assumed, while the number with a causal genetic mutation is smaller under the recessive model assumption.
In summary, if the study population is genetically homogeneous, an independent subjects design provides unbiased genetic estimates, regardless of the modelling approach, and offers the most powerful tests as well. The models that use repeated measurements are even more powerful than those using one single outcome per subject, however. When the study population is stratified, using OLR or OLMM can result in spurious associations; therefore, in order to control for potential population admixture or stratification, familybased designs are strongly recommended. ALMM is more powerful than QTDT_{M }and FBAT statistics.
Discussion
The ALMM method combines the characteristics of both longitudinal data analysis and caseparent design. While repeated measurements of quantitative phenotypes allow for the assessment of the effect of timedependent environmental exposures, the use of the caseparent design with analysis based on ALMM eliminates the potential bias in the estimated coefficients associated with population admixture or stratification, provided that the linear model is correct. Unbiased tests also require valid estimates of standard errors, however [8]. Use of the ALMM model, which allows intercepts and environment effects to depend upon parental mating type, can help to address that issue [14]. Our results show that ALMM represents a valuable methodology for correctly assessing main and geneenvironment interaction effects for quantitative traits in stratified populations. In addition, by taking advantage of the structure of the ordinary linear mixedeffects models, covariates may be included and balanced repeated measurements are not required. This model can be implemented using standard statistical software, including a linear mixed models module.
Since ALMM is based on a caseparent design, ethnicity bias is avoided because all possible genotypes are taken into account, even those that were not transmitted to the affected offspring. Including an indicator variable for the mating type allows one to use different intercepts; thus, differences within and across matingtypes are considered in the genetic effects estimation. In order to account for a potential correlation between the exposure and the mating type, different levels of exposure, depending on the mating type, have been modelled in ALMM through the inclusion of the interaction between such an indicator variable and the exposure. In this manner, situations where, for example, the environmental exposure may depend on the mating type can be assessed without bias. Therefore, population stratification and admixture are no longer sources of estimation bias.
It can be the case that the study population is mixed, although the trait of interest does not vary within the subpopulations. In those situations, ethnicity is not a confounder; thus, genetic effects may be estimated without bias through the use of ordinary regression models. When there is an admixed population, and the exposure of interest does not depend on substructure, the indicator variable for the mating type (α_{1m}) can be omitted in model (13). In the case where Z is randomised, it is known that Z is independent of exposure and genetic background; allowing α_{6m }to depend on m ensures that the effect of treatmentgenetic interaction is not confounded by different responses to treatment in the different subgroups. If it can be assumed that both coefficients do not depend upon m, a simpler model can be fit, which has two clear advantages. First, there are more degrees of freedom, and this is important in small studies, especially if some strata of mating types have few observations. Secondly, if parents are missing, instead of computing E(X) conditional on parental genotypes, we can replace it by E(X) conditional on the sufficient statistics for parental genotype [8].
A disadvantage of caseparent designs is that parental genotypes are not easily accessible for lateonset disorders. In those cases, other familybased designs suggest using siblings rather than parents, although larger sample sizes are required in order to achieve comparable power [21].
It is important to note that, because OLR and OLMM provide biased estimators under population stratification, power comparison against unbiased methods may not be completely fair. Power will be underestimated when the parameter is incorrectly estimated with values that are close to zero, although when the reverse occurs, the power will be magnified. For that reason, it was decided to exclude those methods in power comparisons under population stratification.
For testing both main genetic effect and geneenvironment interaction effects, regardless of the composition of the population, ALMM was found to be more powerful than the twostep modelling approach where AQTDT_{M }and FBAT  or QBATI, in the geneenvironment interaction assessment  were used in the second stage. This is because, while the longitudinal analysis approach takes advantage of both repeated values across time and measurements across people, the twostep procedure does not account for the relative degree of within and betweensubject variability. Nevertheless, there are weighting procedures that account for both sources of variability; thus, the summary statistic obtained in the first step can be adjusted for [22]. Although, methodologically, the linear mixed models represent an adequate approach for longitudinal data analysis, one should not forget about the twostep modelling approach because it represents an intuitively simpler procedure and the opportunity to use existing genetic software. It should be noted that the difference in power between both modelling approaches may strongly depend on the number of repeated measurements, the underlying true effect sizes and the frequency of the missing phenotypes.
It is evident that, with a homogeneous population, OLMM is the most powerful tool. The decrement in the power of ALMM relative to OLMM is related to the lessening of the variability in the genotype due to the centring procedure and to the extra parameters, for each mating type, to be estimated. When the population is not homogeneous, however, these factors should not represent a disadvantage when contrasted with the added advantage of an unbiased estimation.
In summary, in addition to comparing the longitudinal approach against the twostep modelling approach, we have also compared designs using independent subjects against familybased approaches under homogeneous and stratified populations. Assuming no population stratification, ordinary regression methods are valid and more powerful than the other methods. Nevertheless, the familybased approach is strongly recommended when the homogeneous ethnicity in the population is in doubt, in order to achieve unbiased estimators. ALMM now represents a powerful tool for assessing the main genetic effect and geneenvironment interactions on timedependent phenotypes under population stratification.
Appendix 1
R code for ALMM
library(nlme)
mmodel < lme(fef2575~ o3*tx*as.factor(mating_type) + gstm1*tx*o3, random = ~ 1 + o3id, method = "ML", data = base)
where:
fef2575 = outcome
o3 = ozone exposure (timedependent)
tx = supplementation treatment (fixed)
gstm1 = genotype (takes values 0, 1, 2 if additive disease model; or 0, 1 if recessive)
mating_type = vector with values 1, 2, 3, 4, 5, 6 classifying the different mating types.
id = individual identification
R code for the second stage in the twostep modelling approach
QTDT_{M}
qtdtm < lm(slope ~ tx*as.factor(mating_type) + tx*gstm1, data = base1)
FBAT
library (pbatR)
pbat.m(slope ~ tx  gstm1, ped = ped, phe = phe, fbat = "gee", min.info = 10, max.pheno = 1, scan. genetic = "additive")
QBATI
library (pbatR)
pbat.m(slope ~ mi(tx)  gstm1, ped = ped, phe = phe, fbat = "gee", min.info = 10, max.pheno = 1, scan.genetic = "additive")
where:
gstm1, tx and mating_type are defined as above
More details about the PBAT commands can be found in http://biosun1.harvard.edu/~clange/pbat.htm webcite
General FBAT statistic
For N nuclear families, one offspring in the family i and no covariates
where
and E(X_{i}g_{im}, g_{if}) and Var(X_{i}g_{im}, g_{if}) are calculated under the null hypothesis of Mendel's law. That is:
and
where g on the right hand side of these expectations indexes the possible offspring genotypes and P(g) is the probability of a particular genotype given the parents' genotypes, calculated under the null hypothesis. Thus,
If both parents are homozygous, X_{i }= E(X_{i}g_{im}, g_{if}) and Var(X_{i}g_{im}, g_{if}) = 0. Therefore, these triads do not add information to the FBAT statistic and they are referred to as noninformative families.
The test is robust against population stratification, as a result of centring X by its expected value conditional on parental genotypes (g_{im}, g_{if}) assuming Mendel's laws.
The statement that case selection was not based on their genotype information is the only assumption about the ascertainment process.
Since in U, E(Y_{i}) is calculated under the null hypothesis, it can be estimated by . Note that the test statistic is based on the relative size of U with respect to its standard deviation but not on the size of β_{1 }explicitly. Thus, the genetic effect is not directly estimated.
QBATI
This statistic[17] is based on the following regression model:
where:
S_{i }= (g_{im}, g_{if}) Sufficient statistic (parental genotypes) β_{0}(Z_{i}, S_{i}) encodes the dependence of the outcome on the environmental exposure and the parental genotypes
β_{1 }= main genetic effect
β_{2 }= geneenvironment interaction effect
And X_{i }and z_{i }are as defined previously.
Note that there is no coefficient for the environmental effect, as this is subsumed in the intercept β_{0}. Assuming that the environmental exposure is independent of the candidate gene, and conditional on S_{i}, estimators for both β_{1 }and β_{2 }are obtained through the equation:
where:
Under weak regularity conditions, the solution to this equation leads to consistent estimators for β = (β_{1}, β_{2}) which are robust for population stratification.
The test statistic for the geneenvironment interaction has the same form as the original FBAT statistic given in (12); that is:
where:
with:
and
is an estimate for the main genetic effect under the null hypothesis of no geneenvironment interaction and is a weighted average of the environmental exposures that ensures .
Note that, in (14), the point of attention is on the genetic effect through the main and the geneenvironment interaction. In other words, the parental genotype and the environment main effect are not of direct interest for estimation. In this sense, the test of H_{0}: β_{3}= 0 based on model (9) may be thought of as an equivalent test to QBATI
Table S3. Table S3a. Bias results for main genetic effect assessment comparing ordinary statistical methods (OLR and OLMM) with familybased methods (AQTDT_{M }and ALMM) under homogeneous (HP) and stratified (SP) m populations. Each time, n cases were simulated with parameters β_{1 }= α_{5 }= 0.5. Simulations are based on the recessive genetic model. † = number that identifies each model in the paper.
Table S4. Table S4a. Empirical power results for main genetic effect assessment comparing ordinary statistical methods (OLR and OLMM) with familybased methods (AQTDT_{M}, FBAT and ALMM) under homogeneous (HP) and stratified (SP) populations. Each time, n cases were simulated with parameters β_{1 }= α_{5 }= 0.5. Simulations are based on the recessive genetic model. † = number that identifies each model in the paper.
Endnotes
*Usually the variant allele is the less frequent.
†There are six different mating types: AAXAA, AAXAa, AAXaa, AaXAa, AaXaa, and aaXaa. When both parents are homozygous, the observed and conditional expected genotypes are equal and it is said that the family is not informative.
‡Although allele A was not the less frequent, it was considered the variant allele because, in the original study, children with genotype AA were classified as GSTM1 null (no copy), and those with genotypes aA (one copy) and aa (two copies) were considered GSTM1 positive.
Acknowledgements
Dr London was supported by the Division of Intramural Research, National Institutes of Health, Department of Health and Human Services (ES049019). Dr Laird was supported by the National Institute of Mental Health. The authors thank Douglas Dockery, ScD and Diane Gold, ScD from the Environmental Health Department, Harvard School of Public Health, for their invaluable suggestions and comments.
References

Martinez FD: CD14, endotoxin, and asthma risk: Actions and interactions.
Proc Am Thorac Soc 2007, 4:221225. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Hunter D: Geneenvironment interaction in human diseases.
Nat Rev Genet 2005, 6:287298. PubMed Abstract  Publisher Full Text

Li CC: Population subdivision with respect to multiple alleles.
Ann Hum Genet 1969, 33:2329. PubMed Abstract  Publisher Full Text

Deng HW, Chen WM: The power of the transmission disequilibrium test (TDT) with both caseparent and controlparent trios.
Genet Res 2001, 78:289302. PubMed Abstract  Publisher Full Text

Choudhry S, Seibold MA, Borrell LN, Tang H, et al.: Dissecting complex diseases in complex populations: Asthma in Latino Americans.
Proc Am Thorac Soc 2007, 4:226233. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Almasy L: Introduction: Methods for detecting genotype X environment interaction.
Genet Epidemiol 2001, 21:S817S818. PubMed Abstract

Tian C, Gregersen PK, Seldin MF: Accounting for ancestry: Population substructure and genomewide association studies.
Hum Mol Genet 2008, 17(R2):R143R150. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Rabinowitz D, Laird N: A unified approach to adjusting association tests for population admixture with arbitrary pedigree structure and arbitrary missing marker information.
Hum Hered 2000, 50:211223. PubMed Abstract  Publisher Full Text

Gauderman WJ, Macgregor S, Briollais L, Scurrah K, et al.: Longitudinal data analysis in pedigree studies.

Romieu I, SienraMonge JJ, RamírezAguilar M, MorenoMacías HR, et al.: Genetic polymorphism of GSTM1 and antioxidant supplementation influence lung function in relation to ozone exposure in asthmatic children in Mexico City.
Thorax 2004, 59:810. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Laird NM, Horvath S, Xu X: Implementing a unified approach to familybased tests of association.
Genet Epidemiol 2000, 19:S36S42. PubMed Abstract  Publisher Full Text

Allison DB: Transmissiondisequilibrium tests for quantitative traits.
Am J Hum Genet 1997, 60:676690. PubMed Abstract  PubMed Central Full Text

Spielman RS, McGinnis RE, Ewens WJ: Transmission test for linkage disequilibrium: The insulin gene region and insulindependent diabetes mellitus (IDDM).
Am J Hum Genet 1993, 52:506516. PubMed Abstract  PubMed Central Full Text

Ewens WJ, Li M, Spielman RS: A review of familybased tests for linkage disequilibrium between a quantitative trait and a genetic marker.

Gauderman WJ: Candidate gene association analysis for a quantitative trait, using parentoffspring trios.
Genet Epidemiol 2003, 25:327338. PubMed Abstract  Publisher Full Text

Laird NM, Lange C: Familybased designs in the age of largescale geneassociation studies.
Nat Rev Genet 2006, 7:385394. PubMed Abstract  Publisher Full Text

Vansteelandt S, Demeo DL, LaskySu J, Smoller JW, et al.: Testing and estimating geneenvironment interactions in familybased association studies.
Biometrics 2008, 64:458467. PubMed Abstract  Publisher Full Text

Diggle PJ, Liang KY, Zeger SL: Analysis of longitudinal data. Oxford University Press, Inc., New York, NY; 1994.

London SJ, Romieu I: Gene by environment interaction in asthma.
Annu Rev Public Health 2009, 30:5580. PubMed Abstract  Publisher Full Text

Fitzmaurice G, Laird NM, Ware J: Applied Longitudinal Analysis. John Willey and Sons, Inc., Haboken, NJ; 2004.

Abecasis GR, Cardon LR, Cookson WO: A general test of association for quantitative traits in nuclear families.
Am J Hum Genet 2000, 66:279292. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Korn EL, Whittemore AS: Methods for analyzing panel studies of acute health effects of air pollution.
Biometrics 1979, 35:795802. PubMed Abstract  Publisher Full Text