For example, for a continuous outcome, the link function g(x) = x is equivalent to standard linear regression, and for a binary outcome, the link g(x) = log(x/(1 – x)) is equivalent to logistic regression. For ease of exposition, the genotype is assumed to be coded additively (0, 1, or 2 copies of the minor allele) and the exposure is assumed to be binary, so that is a scalar. The results are easily extended to categorical exposure with three or more levels, for example.

We note that for continuous exposures, both standard and joint tests can be invalid when the exposure main effect is mis-specified [8, 9]. Using the Huber-White robust ‘sandwich’ estimate of the variance-covariance matrix of yields valid tests even when the continuous exposure is not modeled accurately. These robust covariance matrices are currently available from some software packages (e.g. ProbABEL [10]), but not others (e.g. PLINK [11]). Alternatively, breaking the continuous variable into categories (e.g. a binary indicator for high versus low exposure) also yields valid tests [8].

Following the exposition in van Houwelingen et al. [12] and assuming the sample size in every study is large enough so that is multivariate normal with variance-covariance matrix Σ_i, we can write the log likelihood for the observed as:

One can solve for the maximum likelihood estimate using the score estimating equation

This leads to the weighted least square solution:

Given , the Wald test of the joint null β = 0 is , where is the Fisher information.

The score test also has a simple form:

Both of these tests have a non-central χ² distribution with 2 d.f. under the null hypothesis. (The assumption of normality is used here to motivate the estimation procedure; the test remains valid if the estimates and Σ_i are unbiased but not normally distributed. Due to the consistency and asymptotic normality of maximum likelihood estimates, the study-specific estimates are likely unbiased and normally distributed when the individual sample sizes are ‘large enough’. Care should be taken when some studies have relatively small sample sizes, e.g. under 100 subjects.)

Although in principle this procedure yields valid and efficient tests of the joint null, in practice it may be difficult to implement, as it requires the covariance of and . Standard packages for the analysis of GWAS data such as PLINK [11] or glu (http://code.google.com/p/glu-genetics/) do not report this covariance without custom modification. (The --robust option in ProbABEL [10] will produce estimates of the variance-covariance matrix.) To avoid this complication, one might conduct a stratified analysis in each study, i.e. estimate the genetic effect in the exposed subjects and the effect in unexposed separately. Since these estimates are uncorrelated (they are constructed using non-overlapping sets of subjects), the score test simplifies to: where and

Note that expression [C] is simply the sum of the usual fixed-effect (inverse-variance weighted) tests for .

To demonstrate this approach, we used genotype data from five GWASs conducted using European-ancestry subjects: case-control studies of breast cancer (BRCA), type 2 diabetes (T2D) and coronary heart disease (CHD) in the Nurses’ Health Study (NHS), and case-control studies of T2D and CHD in the Health Professionals Follow-up Study (HPFS) [13, 14]. For each subject in each GWAS, a continuous phenotype Y was simulated as a function of two single nuclear polymorphisms (SNP) and one binary environmental exposure: where G1 was the (observed) count of minor alleles for rs505922 (frequency from 34.1 to 35.8% in these studies), G2 was the (observed) count of minor alleles for rs1219648 (frequency from 39.6 to 42.0% in these studies), the exposure E was a Bernoulli 0-1 variable with expectation 0.3, and the residual variation was normally distributed with mean 0 and standard deviation σ. E and were independent of each other as well as G1 and G2; since rs505922 is on chromosome 9 and rs1219648 is on chromosome 10, G1 and G2 were also effectively independent. We chose this model to examine the performance of the joint meta-analysis (and the marginal and standard interaction meta-analyses) in two situations: one where the effect does not differ across the exposure strata, and one where it does. The two SNPs were chosen arbitrarily, but are representative of the markers present on genome-wide genotyping platforms (rs505922 is on the Affymetrix Axiom and Illumina HumanHap550 arrays, among others, and rs1219648 is on the Affymetrix 6.0, Illumina HumanHap550, and Illumina OmniExpress arrays, among others).

To illustrate the validity of the meta-analytic joint test, we conducted a GWAS of a single realization of Y. This GWAS included 2,543,290 genotyped or imputed SNPs distributed along the 22 autosomal chromosomes that were available for all of the five studies. Because most of these SNPs are not in linkage disequilibrium with G1 or G2 and hence are not associated with the simulated phenotype, results from tests of the phenotype can be used to estimate the type I error rate of the meta-analytic approaches. For this example, the parameters b = (b_E, b_G1, b_G2, b_G2E)’ and σ were chosen so that the marginal test of G1 would have high power at the 5 × 10⁻⁸ level in the total sample (n = 10,161), and the joint gene-environment test of G2 would have high power in the total sample, while the marginal test of G2 would have modest or low power. The values b = (0.05, 0.12, 0.05, 0.10)andσ = 1.00 satisfied these criteria. Under this model, the genetic variants explained 0.9% of the variation in Y (0.7% for G1 and 0.2% for G2). We estimated subsequently the power of the meta-analytic joint test (and the marginal and standard interaction meta-analyses test) to detect G1 and G2 in this specific situation, by generating 100,000 realizations of Y with these parameters.

To further exemplify the power of the meta-analytic joint test relative to meta-analytic marginal and standard 1 d.f. interaction tests across a broad range of situations, we simulated 1,000 realizations of Y under each of 320 different parameter settings. We kept the values of b_G1, and σ fixed at their values in the previous simulation (0.12 and 1.00, respectively), while we varied b_E, b_G2 and b_G2E across a grid defined by b_E {0.05, 0.1, 0.5, 1, 2}, b_G2 {0.01, 0.04, 0.07, 0.1}, and b_G2E {−0.15 to 0.15 by 0.2 step}.

As estimated by simulation, the power of the marginal test for G1 in this situation was greater than 98.91%, and the power of the joint test for G1 was 97.72%. On the other hand, the power to detect G2 using the joint test was 69.92%, while the observed power of the marginal test was 48.81%. The power of the standard test for G2-E interaction was only 0.86%. As expected, since we did not simulate a non-multiplicative joint effect for G1 and E, none of 100,000 realizations of the standard G1-E interaction test were significant at the 5 × 10⁻⁸ level.

Figure ​Figure22 shows how the power of these tests changes as more studies are included in a meta-analysis. Generally, power increased non-linearly in the number of studies (reflecting the increase in the total sample size), while the power for any individual study was below 1% for all tests of G1 and G2 (from 0.05% for the marginal test of G2 to 0.002% for the standard G2-E interaction test, and from 0.5% for the marginal test of G1 to 0% for the standard G1-E interaction test). In particular, using fewer than all five available studies led to a concerning decrease of the power for the marginal and joint tests of G2. The removal of one study nearly halved power to detect the Y-G2 association using either of these tests. (The standard G2 × E test had very low power even when using all five studies.) Tests of G2 using less than four studies had less than 10% power. The decrease in power to detect G1 when using fewer studies was smaller than for tests of G2: the power of the marginal test for G1 was 85% when using four studies and 45% when using three studies.

Figure ​Figure33 shows the power of the three tests for G2 observed when b_E = 0.05 and varying the values of b_G2 and b_G2E. The power of the joint test to detect G2 is often greater than that of the marginal test. The maximum difference between the power of the joint test and marginal test across these simulated scenarios is 55.46%. The minimum difference, on the other hand, is never more than −7.93%, suggesting that when the joint test has less power than the marginal test its power is not that much smaller. Since the joint test and the standard test of interaction take the environmental effect into account, variation of b_E had no effect on these two tests. However, the power of the marginal test of G1 and G2 decreased for values of b_E over 0.1, and the difference of power of the joint test over the marginal test was increasing in the same proportion. Hence, the joint test of G2 was more powerful than the marginal test of G2 in all situation where b_E was over or equal to 0.5 and the joint test of G1 also was more powerful than the marginal test of G1 when b_E was over 1.0.

Some of the biggest power gains for the joint test come in situations when the main effect of G2 (b_G2) was small and the G2 × E interaction parameter (b_G2E) was large and in the same direction. For example, when b_G2 = 0.04 and b_G2E = 0.13, the power of the joint test was 83%, while the power of the marginal test was 46%. Large gains in power were also observed when the interaction parameter had the opposite sign as the main effect of G2, and consequently the marginal association between G2 and Y was weak. We note that for some of these situations the standard test of interaction has greater power than the marginal test, but it still can have considerably less power than the joint test.

The standard test of gene-environment interaction always showed very low power (under 20% in all simulations), being slightly more powerful than the two other approaches only in some extreme cases, when gene effect was quite low and the interaction effect quite high and opposite to the gene main effect (e.g. b_G2 = 0.04 and b_G2E ≥ −0.13).

Similar to the single-study version of the joint test [5], the meta-analytic test has appropriate type I error rates under the null hypothesis of no association between the tested marker and a phenotype in any stratum defined by the environmental exposure, as long as the marker and the phenotype are uncorrelated in the sampled population [9, 16]. In many situations, including models where the genetic effect in each environmental strata was not null and in the same direction, this meta-analysis of joint test has greater power than either the meta-analysis of standard test for departures from a multiplicative odds ratio model, or the meta-analysis of marginal test for association, making it an attractive test to discover markers which may have been missed by meta-analyses focusing on marginal marker-trait associations (e.g. most meta-analysis of GWAS studies to date).

We have shown that sample sizes far exceeding those of a typical individual GWAS will be needed to reliably detect genes with subtle gene-environment interaction patterns. In the first scenario we simulated, no single study had greater than 1% power to detect a marker whose effect differed across environmental strata (G2). On the other hand, the power of a meta-analysis combining five studies (more than a three-fold increase in sample size relative to the largest individual study) was 58.27 times greater than that of the largest single study; the power of the meta-analysis was greater than 80% for many situations of interest. The large increase in relative power for a smaller increase in relative sample size is due to the non-linear nature of the power curve as a function of sample size. This suggests that meta-analyses incorporating gene-environment interactions may cross a sample-size ‘tipping point’ and identify markers that have gone undetected in individual studies. In the models we simulated, the combined genetic effects of G1 and G2 explained between 0.7 and 1.3% of the total variation in Y, which leads to a high power of the joint and marginal test in the total sample of 10,161 subjects. For a smaller genetic effect, a larger sample size will be required to achieve similar power.

Simulations across a range of genetic main effects and gene-environment interaction effects confirm that in a meta-analysis context, the joint test can be more powerful than either the marginal test of the gene or the standard test of interaction when the genetic effect is weak in one exposure stratum but strong in another, or when the main and interaction effects are in opposite directions (so that the marginal association between the tested marker and the phenotype is weak). However, the latter situation may be of less practical interest since (depending on the context) it may not be biologically plausible; there are few empirical examples of such ‘crossover’ interactions in observational studies [17]. In the range of models we tested, the joint test was often more powerful than the marginal test of the gene, with an absolute gain of power over 50% in the best cases. In situations where the marginal test had better power than the joint test, the difference was small (less than 10% absolute power). The standard test for interaction had low power in all situations we considered; it had slightly greater power than the other two tests only in a few rare, extreme situations. Although we simulated a continuous trait, because the meta-analysis procedure relies on summary parameter estimates (which could be slopes from a linear regression or log odds ratios from a logistic regression), qualitatively similar results should hold for binary disease traits.

Aside from increased power, another potential advantage for the joint test is that it is scale-free. The power and interpretation of tests for the usual gene-environment product term (β_{g e} = 0 in [A]) depend on the link function g(x) [18, 19]. In the context of gene discovery, where the primary goal is establishing associations between markers and a phenotype, no particular scale may have intrinsic interest, and the scale that yields the highest power for the standard interaction test will not be known a priori. For saturated codings of the genetic marker and environmental exposure, the power of the joint test does not depend on the link function for the working model fit to the data [A], but only depends on the true relationship between the marker, exposure and phenotype. However, in situations where departures from additivity on a particular scale are the main interest – for example, when testing for ‘public health interactions’, where the impact of reducing (or increasing) the prevalence of an exposure is larger among genetically susceptible individuals [20] – then the joint test is inappropriate, as it is not a test of interaction per se.

The meta-analysis procedure we propose – [C] – will be useful in typical situations where it is logistically difficult or impossible to share and combine individual-level data. It is easily implemented using standard software for the analysis of GWASs. However, because this particular procedure requires analyses to be conducted within strata of a categorical exposure, in some situations it may lose power relative to analyses that assume a particular form for the interaction term (e.g. ordinal G times continuous E). Moreover, although it does not require that the genetic marker and environmental exposure are independently distributed in the sampled population for validity, it cannot leverage this assumption to increase the power of the individual-study joint tests, as others have proposed [2, 4].