The Association between SNPs and a Quantitative Trait: Power Calculation.

A power calculation for a study with a quantitative outcome requires information on the outcome distribution under the alternative hypothesis. Researchers face challenges when they concisely specify alternative distributions in genetic studies because power depends on genotype frequencies and the average effect of each genotype. In GWAS, investigators evaluate hundreds of thousands of associations; therefore it is unrealistic to specify gene frequencies and gene effects for each test and some simplification is needed. Software packages are available to calculate power, but many of them have limited flexibility and / or may have a steep learning curve. In this review, we describe to researchers and graduate students the essentials of a power calculation for testing for an association between a quantitative trait and genotypes. In addition, we provide them with the codes of the different available software packages-free and commercial-to calculate this power. The calculations can be carried out using virtually any computer language that computes the cumulative distribution function of a non-central F-distribution.

INTRODUCTION

Researchers conduct Genetic-Wide Association Studies (GWAS) to locate genetic variants (usually SNPs: single nucleotide polymorphisms) responsible for common and complex diseases. An insufficient sample will result in a lack of statistical power, which can prevent researchers from identifying SNPs related to such diseases. Therefore, before conducting a GWAS, investigators frequently calculate the ‘required’ sample size to achieve a specified or the statistical power of available samples. Many investigators have no information on how to calculate the required power to conduct GWAS.

In general, power is calculated by specifying a statistical test, e.g., analysis of variance, an effect size to be detected by the test, and the probability of incorrectly rejecting the null (no effect) hypothesis. The researcher wishes to reject the null hypothesis when the alternative hypothesis (effect of specified size) is true. The statistical power of a test is the probability of rejecting the null hypothesis when it is false (accept the alternative hypothesis). In genetic association studies, the power calculation depends on genotype frequencies and the average effect of each genotype. Because GWAS evaluate hundreds of thousands of associations, it is unrealistic to specify gene frequencies and gene effects for each test and some simplification is needed. Herein, we describe how to calculate power to test for an association between a quantitative trait and genotypes by an analysis of variance (ANOVA). To describe the effect size of SNP’s genotypes on quantitative traits, non-centrality parameter (NCP) is used. The NCP defines the effect size in an ANOVA (Lehmann, 1986). Once a NCP is computed, the power calculation is a function call that computes the cumulative probability of a non-central F distribution, a widely available function. Thus, 1) it is easy to compute power for an F-test of phenotype differences among genotypes. 2) It is easy to implement this calculation in many statistical software packages. If the function is unavailable, there are good numerical approximations that are not difficult to program (Abramowitz and Stegun, 1965; Johnson et al., 1995).

In this paper, we explain how to compute a NCP based on genetic concepts of phenotype heritability, additive allele effects, or dominant/recessive allele effects. We illustrate the basic power calculations coded in Mathematica, Matlab, R, SAS and Stata. We outline algorithms to plot power as a function of heritability and implement the algorithm in SAS. We outline algorithms to plot the effect size detected with power of 0.8 as a function of minor allele frequency and implement the algorithm in Matlab. Our approach relates the statistically relevant NCP to genetically relevant parameters to provide better insight into the power of a study to identify genetic variants associated with quantitative traits.

Although there are several stand-alone packages that can compute power (Dupont and Plummer, 1990, 1998; Faul et al., 2007; Feng et al., 2011; Gauderman, 2002; Purcell et al., 2003; Spencer et al., 2009), learning how to use them takes time and they tend to lack flexibility. Our goal is to make the power calculation straightforward for someone who is familiar with any of the common statistical packages.

METHODS

We assume that associations between a continuous phenotype, e.g., blood pressure, and many genetic markers will be tested one at a time by analysis of variance (ANOVA). Table 1 is a generic ANOVA summary of the analysis for a single marker. Several terms of interest to the discussion of the statistical power of these analyses are symbolized in this table. The table assumes that phenotypes and genotypes are measured in n subjects. The degrees of freedom for genotypes, d, will depend on how we evaluate the genotypes, but typically d will be 1 or 2.

Table 1.

ANOVA table testing for changes in phenotype associated with the genotypes of a SNP

Source of Variation	Degrees of Freedom	Sum of Squares	Expected Mean Squared Error	F-ratio
Covariates	p − 1	*	*	*
Genotypes	d	SSgenotypes	σ2 + η2	F = v SS_genotypesd SS_residual
Error	v = n − p − d	SSresidual	σ2
Total	n − 1

The p covariates are incorporated to adjust for non-genetic factors that may affect the phenotype. In the case of blood pressure for example, these might be covariates that adjust for age, blood pressure medication, and population stratification. The same covariate adjustment is applied to the analysis of each genetic marker. So operationally the phenotype measurements are initially adjusted by these covariates and the marker-specific computations are usually run on the residuals after the covariate adjustment. An ANOVA test of ‘no gene association’ tests the null hypothesis that there are no phenotype differences among genotypes, i.e., η2 = 0. Here, η2 symbolizes the effect size, which is a summary of the magnitude of the phenotype differences among genotypes. The power of this F-test can be computed using functions that are widely available in statistical software packages. Power depends on the magnitude of a so-called non-centrality parameter (NCP), λ = η2 ∕ σ2, the degrees of freedom, d and ν, and the significance level, a. In particular, the power of this test is computed as the probability that the F-ratio will exceed a value, Power = Pr[F(d, ν, λ) ≥ c], where F(d, ν, λ) denotes a random variable with a non-central F distribution and c denotes the value of an F-ratio that is statistically significant at the a–level, i.e., c such that a = Pr[F(d, ν, λ = 0) ≥ c]. These two values, Power and c, can be calculated in many software packages, e.g., SAS or Matlab. Thus, it is straightforward to compute power for an F-test of phenotype differences among genotypes when we specify a NCP. Unfortunately a NCP does not provide an insightful interpretation of the effect size in the context of a GWAS that seeks to identify differences in a continuous phenotype associated with single nucleotide polymorphisms (SNP).

NCP as a Function of Heritability

Our initial aim is to relate the non-centrality parameter to the proportion of variation that is explained by the SNP, R2. When the computations are run on the residuals after covariate adjustment, R2 is simply the R-squared associated with adding genotype. With respect to entries in Table 1, this is the observed proportion of the phenotype variation that is explained by the SNP genotypes after adjusting for the covariates,

The expected value of a non-central F-distribution is:

The expected F is a one-to-one function of the non-centrality parameter,λ, and consequently specifyingμis equivalent to specifying the non-centrality parameter, specifically

The ‘true’ proportion of the phenotype variation (after adjustment for covariates) that is explained by the genotypes is denoted π. Replacing the F-ratio in Table 1 with the expectations, yields a relationship that specifies the non-centrality parameter as a function of the expected proportion of the phenotype variation that is explained by the genotypes, i.e.,

Genetic Models

The entries in Table 1 assume a simple model for phenotype heritability where each genotype yields a phenotype that has a normal distribution with genotype-specific mean and a variance that is the same across genotypes. For our purpose here, it is convenient to assume that there are two alleles, one of which has an allele frequency less than or equal to the other, a.k.a. minor allele. The genotype frequencies of SNPs are usually unknown at the experimental design stage. However for power calculations, we approximate genotype frequencies in terms of their minor allele frequency (MAF), and further assume that genotypes are in Hardy-Weinberg equilibrium. That is, genotype frequencies for the genotype with 0, 1, or 2 minor alleles, {f0, f1, f2}, are approximated by {(1 − q)2,2q(1 − q), q2} where q is the MAF. Under this model, , where μdesignates the mean value of the phenotype distribution for a genotype with ‘a’ minor alleles. In this general setting, the task of specifying μ for every marker is too tedious to be useful, and we recommend the preceding approach for obtaining the NCP. However, there are two genetic models which provide straightforward interpretations.

NCP in an Additive Genetic Model

An additive genetic model assumes that alleles have an additive effect on the phenotype and our assumptions allow us to relate the NCP to the trend in the phenotype means, β, with the number of minor alleles. To model this, we count the number of a specified allele in the genotype, X = 0, 1, 2, and write the phenotype means as μ = μ0 + βx. Note that this formulation implies that μ = μ0 + βx and μ − μ = β(x − x). The test for a trend in the phenotype with increasing count of minor alleles is a test that β = 0. In the context of Table 1, this trend test corresponds to d = 1 and where f is the genotype frequency and β is the slope of the trend (change in phenotype mean per minor allele). These assumptions allow us to relate the effect size, β/σ, to the NCP,

NCP in dominant/recessive models

We can co-opt the previous results to apply to dominant genetic models. If the minor allele is dominant, let and if the minor allele is recessive let with either off these changes to the definition of x, the NCP can be computed using the same formulae as the additive genetic model. However the interpretation of the slope parameter, β, changes to the difference between the phenotype means of the dominant genotypes and the recessive genotype, respectively. Note that the power is the same whether the minor allele is dominant or recessive since β computed under the dominant assumption simply becomes −β under the recessive assumption; both yield the same NCP power is the same whether the minor allele is dominant or recessive since β computed under the dominant assumption simply becomes −β under the recessive assumption; both yield the same NCP.

RESULTS

Power is the probability of a significant ANOVA under the alternative hypothesis, i.e., Pr(F ≥ c) where F denotes a random variable with a non-central F-distribution, F(d, n-p-d, NCP), and c is the value that is significant at the a-level : Critical Value. Power and c can be calculated in many software packages. Figure 1 illustrates the fundamental power calculation coded in several well-known computer packages. The codes in Figure 2 outlines the basic power calculation. The user specifies the type 1 error (significance level), sample size (n), number of adjusting covariates (p), the numerator degrees of freedom (d), and the NCP. Typically a marker SNP has two alleles and a minor allele frequency (MAF) exceeding 5%. In a large sample, n = 1000, at least a few individuals of the three possible genotypes should be present in the sample giving the degrees of freedom for genotypes, d= 2. The code in Figure 1 can be embedded in a loop that varies the NCP as a function of the expected proportion of the variation that is explained by the genotypes (Eqn 1) or as a function of β/σ (Eqn 2). Figure 2 plots power curves based on the expected proportion of the variation that is explained by the genotypes. This figure iterates the preceding code for sample sizes, n = 500, 1000, 2000, 4000, 8000, and NCP defined by Eqn 1 with π = 0 to 0.01 by 0.001. Figure 2 was produced using Matlab as the authors are more adept at plotting using this programming language. SAS code that will reproduce the essential features of Figure 2 are included in the appendix. Because the NCP in Eqn 2 also depends on the MAF, we report the effect size, β∕σ, that has a power of 0.8 as a function of the MAF. Conceptually, this could be computed by adding another loop that varies the MAF and saving the effect sizes with a power of 0.8. Figure 3 plots the effect size detectable with a power of 0.8 as a function of the MAF for the sample sizes used in Figure 2. This was computed using a Matlab program. The calculations can be finessed to express power in the context of a disease locus that is observed via linked markers. The appended Matlab program incorporates an adjustment for the correlation between a disease locus and the markers.

Figure 1.

Computer code in SAS, Matlab, R, Stata, and Mathematica that compute power for an ANOVA (Table 1)

Figure 2.

Power of an ANOVA for several sample sizes as a function of the proportion of the phenotype variation that is explained by the genotypes

Figure 3.

Effect size detectable with a power of 0.8 in the dominant/recessive genetic model as a function of the minor allele frequency for several sample sizes

DISCUSSION

In this paper, we described computer codes that calculate the required statistical power for detecting a significant association between a genetic variant and a quantitative trait in GWAS under additive, dominant and recessive models for the phenotype — Power Calculation for SNPs and quantitative trait Association (PCSQT). These codes, i.e. PCSQT, were implemented in five widely used commercial statistical packages. The other statistical methods do not explicitly include this parameter while calculating study power. Our approach relates the statistically relevant NCP to genetically relevant parameters to provide better insight into the identification of genetic variants associated with quantitative traits. If the NCP is specified it is easy to compute statistical power for an F-test of phenotype differences among genotypes that can easily be implemented in many statistical software packages. The calculations can answer the following question: if a Single Nuclear Polymorphism (SNP) explains % of the variation i.e. Heritability, with sample size n, what is the statistical power of the proposed study?

These codes permit investigators to employ heritability as an effect size, instead of the mean of a quantitative trait, which is often unavailable in genetic studies. To calculate power, we use the following parameters: (2) Total sample size; (2) Heritability (the range is 0~1); (3) type 1 error rate: (4) number of SNPs in the genetic study: (5) number of covariates; 6) linkage Disequilibrium (LD) r2 (Pritchard and Przeworski, 2001): investigators input a LD value between a genetic marker and a hypothesized causative variant. When the LD between the genetic and causative markers is r2, the sample size (N) is increased to be N/r2. In other words, to achieve approximately the same power with the genetic marker as is achieved with the causative variant, the sample size must be increased by a factor of 1/r2. (Pritchard and Przeworski, 2001).

The output information include: (1) statistical power; (2) a family of power curves plots with different heritability and sample size combinations, as shown in Figure 1. The suggested statistical approach has some limitations. First, population stratification is not considered while calculating statistical power using these codes. Nevertheless, it can be added as one of the covariates while estimating the required power calculation. Second, the effect of the interaction between genetic variants, and environmental factors on the power determination cannot be included in the power calculation. In a GWAS that examines a million SNPs, a Bonferroni adjustment for multiple comparisons implies that the critical value, c, is computed for a significance level on the order of 10−8 (Klein, 2007). The numerical algorithm implemented in the software to approximate c may or may not be adequate for the task. For the five packages. the software packages in Figure 1, we could not find documentation for the algorithm or the error in the numerical approximation. Mathematica allows one to specify the number of significant digits; we set this to 20 and use this as a standard in Table 2. Table 2 presents an example of critical values computed to 6 decimal places for examined software packages and where the values are computed, all agree except for α = 10−12. SAS was the only software package that failed to compute the critical values for α > 10−5. Differences among critical values in Table 2 have little practical importance. In practice, the F-distribution will only approximate the actual distribution of the ratio statistic in an analysis of variance (Table 1). The approximation is satisfactory at conventional significance levels such as 0.05 or 0.01, and p-values computed from the F-distribution are generally accepted as ‘good’ by the scientific community. However, very small p-values, such as 10−8, imply a degree of approximation to the distribution of the test statistic that is unlikely to be achieved.

Table 2.

Software packages and Critical value for F-distribution with 1, 1000 degrees of freedom

Significance level	Critical value for F-distribution with 1, 1000 degrees of freedom
	Mathematica	SAS	Matlab	R	Stata
1E-03	10.891865559	10.891866	10.891866	10.891866	10.891866
1E-04	15.259521389	15.259521	15.259521	15.259521	15.259521
1E-05	19.712947049	19.712947	19.712947	19.712947	19.712948
1E-06	24.228934152	.	24.228934	24.228934	24.228933
1E-07	28.794927827	.	28.794928	28.794928	28.794928
1E-08	33.403406313	.	33.403406	33.403406	33.403408
1E-09	38.049531722	.	38.049532	38.049532	38.049530
1E-10	42.730026564	.	42.730026	42.730026	42.730026
1E-11	47.442581611	.	47.442581	47.442581	47.442581
1E-12	52.185519872	.	52.185566	52.185566	52.185520