Due to the rapid advances in genotyping technology, genome-wide association studies with hundreds of thousands of markers are now possible. These large scale genetic studies offer great promise to expedite the discovery of the common genetic variants affecting common diseases [1]. A first step in the analyses is to understand the properties of the massive data sets. Among the most fundamental properties are the proportion of markers without true effects (p₀) and the effect sizes (Δ) of markers with effects. Knowledge of these parameters provides information about how relevant the genotyped markers are for the disease outcome. In addition, these parameters play a role in a variety of applications. For example, estimates of p₀ are commonly used in methods for controlling false discoveries, which is important to prevent spending time and resources on leads that will eventually prove irrelevant. Another example is that knowledge of effect size Δ is important to design follow-up studies that have adequate power to replicate previous findings.

Multiple methods have been proposed to estimate p₀ [2–4]. These estimators tend to make general assumptions about the distribution of the test statistic under the alternative hypothesis, partly because many of them have been developed in the context of microarray research where specific assumptions may be problematic. However, in genetics good approximations for the statistical test statistic distribution are often available. This information can be used to obtain more precise estimators in these applications. Once (a set of) markers have been identified as being associated with the disease, the next objective of interest is typically to estimate the effect sizes. Currently, the most commonly used approach simply estimates the effect size of the significant markers in the same sample that has been used for testing. Due to the effects of sampling error and the presence of false positives, this approach overestimates the effect sizes considerably [5, 6]. Several methods have been proposed to obtain more unbiased effect size estimates such as a simple split-sample method, cross-validation, and bootstrap resampling [7]. However, with all these methods that use the same sample to first declare significance and then estimate the effects sizes of the significant finding, it will remain difficult to obtain estimates that are both precise and unbiased. We therefore proposed a set of related methods (unpublished data) that estimate the average effect size Δ of the 1 - p₀ markers with effects. Because our estimates are not confined to only those markers that are declared significant, they do not suffer from the upward bias caused by sampling fluctuations producing large test statistics in this specific sample. Our methods do not assign effect size estimates to individual markers but estimate the average of all markers with effects. This does not hamper the design of subsequent replication studies. Thus, for any critical value chosen to declare significance, we can calculate the number of markers with effects plus their average effect size among the significant markers to design replication studies.

Our methods include a maximum likelihood (ML) and a quasi-ML (QML) approach that incorporates the test statistic distribution and estimates Δ simultaneously with p₀. In addition, we propose a conservative estimators of p₀ (CON) and a variation of this conservative estimator that adaptively estimates a fine tuning parameter (ADA CON). Neither CON nor ADA CON rely on the test statistic distribution under the alternative but take advantage of the specific knowledge that in large-scale genetic studies the p₀ must be very close to 1. Because these conservative estimators do not consider the distribution under the alternative hypothesis, they cannot estimate the average effect size directly. However, we can still use the point estimate of these conservative methods and include it in a second step in a ML method to estimate the average effect size Δ for our conservative estimators of p₀. We apply our methods to the simulated rheumatoid arthritis (RA) case-control data with 10 k single-nucleotide polymorphisms (SNPs) in Genetic Analysis Workshop 15 (GAW15). We chose a case-control design with SNPs because this is one of the most important designs for mapping the genetic determinants of complex human diseases through genome-wide association studies. We also compared our estimators with four existing estimators. We found two studies comparing multiple and non-overlapping sets of estimators [2, 3]. In these studies, the lowest slope (LOW S) and location based estimator (LBE) showed the most favorable properties and were therefore included here. In addition, estimators developed by Storey (STO) [4] and Storey-Tibshirani (STO-TIB) [8] were included because they may be among the more commonly used estimators.

The maximum likelihood methods (ML and QML) and the conservative methods (CON and ADA CON) are briefly described below. All 100 replicates in GAW15 Problem 3 were analyzed. To create a case-control data set we selected the first sib from family-based data sets as independent cases (N = 1500), and used all individuals in control data sets.

Analyses were done with knowledge of the "answers" of causal markers locations.

Based on the p₀ estimate in our new methods, the average number of total causal markers with main effects was 18. The average numbers of causal markers in the LOW S, STO-TIB, STO, and LBE were 21, 150, 256, and 0, respectively. Clearly, STO-TIB and STO overestimated the number of effects and LBE underestimated the number of effects. It is also important to note that standard errors of the estimates were about 100 times larger for STO-TIB and STO, implying that the number of markers was estimated very imprecisely.

The second part of Table 1 shows results for the estimated average effect size Δ. The ML methods estimate Δ and p₀ simultaneously. The other estimators do not consider the distribution under the alternative hypothesis, and can therefore not estimate the average effect size directly. However, in these cases we can still use the point estimate $p_0^{\ast }$ obtained with these estimators and include that in the a maximum likelihood method that finds $\stackrel{⌢}{\Delta }$ by maximizing ℓ(m - m$p_0^{\ast }$, Δ). In cases where the point estimate $p_0^{\ast }$ equals 1, the effect size cannot be estimated. In these scenarios Δ was treated as "missing". Results showed that the estimated average effect size was 0.083 in all four new methods in which the ML method was slightly more precise. The estimated average effect size was less precise and considerably lower with STO and STO-TIB, reflecting the downward bias and larger standard deviation in these p₀ estimates.

Results illustrated that all of our four new estimators have favorable properties in terms of the standard deviation with which p₀ is estimated. The ML and QML estimators have the additional advantage that they provide a direct estimate of average effect size Δ. Because the point estimates of p₀ in both CON and ADA CON methods are very similar to that in the ML and QML methods, the average effect size is expected to be similar across methods. This is important because these two parameters are somewhat intertwined and the estimate of the average effect size helps the interpretation of the p₀ estimate. For example, without this effect size estimate, it is unclear whether the estimated numbers of causal markers have very small or large effects.

On average, the ML method performed slightly better than the QML method. Furthermore, we found in other simulations that the QML estimator can be unstable. In general, the ML method may therefore be the method of choice. Results also showed that the CON method performed well and was even slightly more precise than the ML estimators. One reason is that the CON method only estimates a single parameter, whereas the ML methods estimate two parameters. However, this observation is also consistent with previous simulations showing that in less optimal conditions (small sample sizes and small number of markers), the CON method can be more robust. Indeed, as another example of its relative robustness, the CON method performed equally well when the four markers with extremely large effects were included but the ML estimators became somewhat less precise.

Linkage disequilibrium causes test statistics between markers to be correlated. Extensive simulations were performed to examine the impact of such correlated tests on our estimates of the p₀ and Δ (data not shown). Results demonstrated that correlated tests mainly increase the variance of these estimates but did not introduce bias. This makes intuitive sense and essentially mimics other scenarios where certain statistics (e.g., mean) are estimated with correlated observations.

Further improvements and extensions of the proposed methods are conceivable. An example involves work we are currently doing to estimate the distribution of effect sizes. The extension essentially consists of conditioning on the number of markers with effects and then maximizing the likelihood L(Δ|m₁). Thus, we start with estimating the largest effect in the data set, then the second largest, continuing until the estimated effect sizes become (very) small. Another example is that in case-control studies, population stratification can cause spurious associations between marker alleles and disease status when both disease prevalence and allele frequencies differ among subgroups. Using the principle of genomic control [1, 13, 14], our estimators can be further adapted to obtain estimates of p₀ and Δ that take stratification into account.