Power estimation and sample size determination for replication studies of genome-wide association studies

Background Replication study is a commonly used verification method to filter out false positives in genome-wide association studies (GWAS). If an association can be confirmed in a replication study, it will have a high confidence to be true positive. To design a replication study, traditional approaches calculate power by treating replication study as another independent primary study. These approaches do not use the information given by primary study. Besides, they need to specify a minimum detectable effect size, which may be subjective. One may think to replace the minimum effect size with the observed effect sizes in the power calculation. However, this approach will make the designed replication study underpowered since we are only interested in the positive associations from the primary study and the problem of the “winner’s curse” will occur. Results An Empirical Bayes (EB) based method is proposed to estimate the power of replication study for each association. The corresponding credible interval is estimated in the proposed approach. Simulation experiments show that our method is better than other plug-in based estimators in terms of overcoming the winner’s curse and providing higher estimation accuracy. The coverage probability of given credible interval is well-calibrated in the simulation experiments. Weighted average method is used to estimate the average power of all underlying true associations. This is used to determine the sample size of replication study. Sample sizes are estimated on 6 diseases from Wellcome Trust Case Control Consortium (WTCCC) using our method. They are higher than sample sizes estimated by plugging observed effect sizes in power calculation. Conclusions Our new method can objectively determine replication study’s sample size by using information extracted from primary study. Also the winner’s curse is alleviated. Thus, it is a better choice when designing replication studies of GWAS. The R-package is available at: http://bioinformatics.ust.hk/RPower.html.

Hence, all findings from GWAS need to be verified. Replication study is a commonly used approach to verifying positive findings [3,4]. If an association between one specific SNP and a certain disease has been identified in the primary study and confirmed in the replication study, we usually treat this association as true positive with a high confidence. If an association identified in the primary study cannot be confirmed in the replication study, we often suspect that it is a false positive.
The power of replication study is crucial in this validation process. If the replication study is underpowered, then the positive findings will have a low chance to be replicated. It's essential to design a replication study with enough statistical power.
How to estimate the power of a replication study in the design phase?
Traditionally, a replication study is regarded as another independent primary study. Thus, the same power calculation in the original primary study is used. For the associations identified in the primary study, a minimum effect size needs to be specified. Then, the underlying alternative distribution of test statistics is assumed to have specified effect size. The major limitation of this traditional power calculation method is that the specification of the effect size is subjective and may cause bias. Besides, no information from primary study has really been used.
One may think to plug the observed effect sizes from the primary study in the power calculation of the replication study. This power estimation approach doesn't need to specify any parameters. Since only significant associations are considered in the replication study, the observed effect sizes for those associations will tend to be overestimated [5]. This phenomenon is known as the "winner's curse" [6], which makes the estimated powers tend to have higher values.
A lot of methods have been proposed to overcome the winner's curse in effect size estimation. An incomplete list includes conditional maximum likelihood estimation (CMLE, [7][8][9]), bootstrap [10], full Bayesian method [11] and Empirical Bayes method (EB, [12]). Since power function is usually not a linear function of effect size, the estimators obtained by simply plugging those bias-corrected effect sizes in power calculation may not achieve the best performance.
Moreover, there are two other challenges in designing replication study: 1. Due to the nonlinear nature and restricted range (limited to [0, 1]) of power function, the distribution of power is usually non-normal when effect size is normally distributed (illustrated in Fig. 1). The interval estimation of the power should consider the non-normality.
2. Since the power values of different associations in the primary study are different, a summary value is needed to determine the sample size of replication study.
This paper aims at addressing the above challenges. Our contributions are listed in the following: 1. For each association identified from the primary study, an EB based method is proposed to estimate its power in the replication study. 2. Due to the non-normality of the estimated power and the inaccuracy of the hyperparameters estimation, a novel interval estimation method combining Monte Carlo sampling and Bootstrap is proposed to estimate the corresponding credible interval of each association's power in the replication study. 3. The average power of the discovered true associations is used for determining the sample size of replication study. An weighted average method is proposed to estimate the average power. Our proposed interval estimation method can also be used to construct the credible interval of the average power. 4. Only the summary statistics of the primary study are needed when using our proposed method to design a replication study. This feature is helpful since summary statistics are more accessible than individual-level genotype data due to the privacy issue and other constraints.
The rest of this paper is organized as follows. In section 'Methods', we will introduce the Bayesian framework to estimate the power of replication studies. We will prove that Bayesian predictive power is immune to the winner's curse. Then we will present how to estimate the power with two-component mixture prior under the Bayesian framework. We will also give the details about estimation of hyperparameters, interval estimation and the estimation of average power. In section 'Results and discussion', we will first use simulation results to demonstrate that our EB based method is better than other plug-in based estimators in terms of overcoming the winner's curse and providing higher estimation accuracy. We will also demonstrate that the coverage probability of given credible interval is well-calibrated. Then we will show the sample sizes determined to replicate findings of 6 diseases from Wellcome Trust Case Control Consortium (WTCCC) [13], which are much higher than the sample sizes estimated by plugging observed effect sizes in the power calculation formula. The increased sample sizes are reasonable due to the winner's curse. In the same section, we will discuss limitations of current modeling and estimation approach. Section 'Conclusions' concludes the paper. When log-odds ratio is normally distributed, the power of replication study possesses non-normality. Assume the effect size follows normal distribution N(0.15, 0.05 2 ). The standard error of log-odds ratio is σ (2) = 0.05 in the replication study. The significance level of the replication study is α (2) = 5 × 10 −3 . The histogram of the power in the replication study possesses non-normality

Methods
We use parenthesized superscript "(j)" to denote primary study (j = 1) and replication study (j = 2). For example, we denote the sample size in the primary study as n (1) . The sample size in the control group and case group are n (1) 0 and n (1) 1 , respectively. The total number of SNPs genotyped in the primary study is m. Among those genotyped SNPs, the proportion of the SNPs having no association with the disease (null SNPs) is π 0 (0 ≤ π 0 ≤ 1).
In both the primary study and the replication study, a contingency table can be created as in Table 1 for each genotyped SNP. With the contingency table, the logarithm of the observed odds ratio reads: ( 1 ) The true value of the log odds ratio μ is usually unknown. The asymptotic standard error of μ (j) can be approximated using Woolf 's method [14], To test whether there is an association between the SNP and the disease, two hypotheses are set up: Wald test can be used to examine whether the null hypothesis should be rejected. The test statistic is z (j) = μ (j) /σ (j) . The significance levels in the primary study and the replication study are fixed to α 1 and α 2 , respectively. Two-sided test is used in primary study. The rejection region is |Z (1) | > z α 1 /2 (We use uppercase letter to indicate a random variable), where z u (0 ≤ u ≤ 0.5) is the upper u quantile of the standard normal distribution N(0, 1). For a replicated association, the test statistics in two studies should be consistent with the same sign. Hence, the test can be regarded as one-sided test in replication study. The rejection region is sgn z (1) where the sign function is ( 4 )

Bayesian predictive power
For an association identified in primary study, the power function in replication study is defined as A traditional power calculation method needs to specify a minimum detectable effect size μ min first. Then, the power of replication study is β (2) (μ min ). Consequently, the power can be used to determine the sample size.
To incorporate information from primary study, the post-hoc method estimates the power of each association by plugging the observed effect size in Eq. (5), i.e. β (2) μ (1) . This approach is widely criticized for the reason of the winner's curse. The estimated power is biased upward since only significant associations are selected in the replication study. To address this problem, a lot of methods have been proposed to overcome the winner's curse in effect size estimation [7][8][9][10][11][12]. Conditional maximum likelihood estimation (CMLE) is the most commonly used type [7][8][9]. In CMLE, the effect size is estimated by maximizing the likelihood conditioning on rejected region, i.e. μ (1) Please notice that, although the selection bias can be reduced using estimator which can adjust estimated effect size, no unbiased estimator exists [11]. With estimated effect size, the power of replication study can be obtained by using β (2) μ (1) CMLE . The plug-in based power estimator is not optimized in terms of minimizing Bayes risk.
The Bayes risk averages the loss function over both sample space and parameter space. In terms of overcoming the winner's curse, the Bayes risk R( θ) should be defined conditioning on rejected region and alternative hypothesis, where θ is the power estimator of replication study. Inspired by the proof in [12], we can show that Bayesian predictive power η (2) [15] is the estimator minimizing R( θ) (please see Appendix for detail). The Bayesian predictive power reads which takes the average of all power function values among all possible μ values given observed z (1) . We will provide a detailed formula of the Bayesian predictive power under one specific prior in the following subsection.

Two-component mixture prior
In each study, the observed log odds ratio μ (j) asymptotically follows normal distribution N μ, σ (j) 2 . The underlying true value of the effect size μ is often unknown. It is widely suspected that a large proportion of SNPs with small effect sizes are associated with complex diseases [16,17]. We use Gaussian prior to depict this pattern of the associated SNPs. For all SNPs, we use the following two-component mixture prior to describe their effect sizes: where δ 0 is the distribution with point mass on zero and σ 2 0 is the variance of the effect sizes in associated SNPs. With this prior, the posterior distribution of effect size μ under H 1 is The Bayesian predictive power of replication study is (Detail in the Appendix): where z * = λ μ (1) /σ (2) , σ * = 1 + λ σ (1) σ (2) 2 and (x) is the cumulative density function (cdf ) of N(0, 1). By substituting observed allele frequencies from the primary study into Woolf 's method, σ (2) can be approximated as There is an unknown hyperparameter σ 2 0 in the calculation of Bayesian predictive power. In the following subsection, we will present how to estimate σ 2 0 with Empirical Bayes approach.

Hyperparameter σ 2 0
In Empirical Bayes' thinking, we can estimate σ 2 0 by taking advantage of the shared structure of the effect size's dis-tribution among all SNPs, which can be seen from Eq. (9). The estimator of σ 2 0 is (see Appendix for detail): There are two extreme cases in the above estimation: 1. If the null hypothesis is valid, then all SNPs follow a standard normal distribution with variance equal to , the sample variance is no bigger than one as in the null hypothesis case, we will have σ 2 0 = 0. In this case, the result of our EB based power estimation method will degenerate to type I error rate, which is the probability that the identified association can be replicated even when the association doesn't exist.

When
In this case, the above shrinkage coefficient will degenerate to λ = 1. The shrinkage effect in our EB based method will disappear.
Noticed that there is another unknown hyperparameter π 0 in calculating σ 2 0 . The estimation of the proportion of true null hypotheses π 0 has been extensively studied [18][19][20]. Here we just choose Storey's method [18] for the simplicity of implementation. Let's denote the number of SNPs with p-value> γ as m + (γ ) in the primary study. Then π 0 can be estimated by using There is a bias-variance tradeoff in tuning γ . An automatic procedure is proposed in [18] without tuning γ : A natural cubic spline will fit to evaluated values with different γ , then π 0 is the spline's value at γ = 1. σ 2 0 can be calculated by plugging π 0 in Eq. (13). By plugging σ 2 0 into Eq. (11), an EB based estimator of the replication study's power can be obtained, which is denoted as η (2) EB . The corresponding credible interval can be constructed, which is presented in the following subsection.

Credible interval
From Eq. (10), the posterior distribution of log odds ratio μ under alternative hypothesis H 1 is a normal distribution. Figure 1 shows the histogram of power values when μ is normally distributed. The shape of the histogram indicates the non-normality of the calculated power. Hence, the asymptotic approach based on normal distribution theory is not appropriate in the interval esti-mation of the replication study's power. The construction of the credible interval should consider the non-normality. We propose to use Monte Carlo sampling to construct the credible interval of β (2) (μ). The credible interval is constructed with known hyperparameters σ 2 0 . Since estimation error will occur in estimating σ 2 0 , the constructed credible interval will have smaller coverage probability than nominal level. To incorporate the variance of the estimator σ 2 0 , a method combining Monte Carlo sampling and Bootstrap is proposed. The test statistics from the primary study z (1) will be resampled N 1 times with replacement. For each run, σ 2 0 is re-estimated. Monte Carlo sampling is used to generate N 2 power values with each re-estimated σ 2 0 . The credible interval is constructed among all N 1 N 2 sampled power values.

Average power
Usually, multiple associations are identified in primary study. To design a replication study, a summarized value reflecting the average power of all associated SNPs is needed. A direct thinking is to average power among the identified SNPs with underlying true associations, which reads β (2) where S is the index set of the associated SNPs identified from primary study and |S| is the cardinality of S. The subscript i means that the quantity is evaluated for SNP i.
Since the index set S is unknown, we propose to use weighted average of the estimated powers η (2) EB . The local true discovery rate (ltdr) of each SNP is the posterior probability of being associated SNP given observed statistics, which is complementary to local false discovery rate [21]. We use ltdr as weight in the estimation. The estimated average power is where T is the index set of the SNPs identified from the primary study. The local true discovery rate of the primary study can be calculated as (See Appendix for detail): where φ(x) is the probability density function (pdf ) of  N(0, 1).
By setting the estimated average power larger than a threshold, e.g.η (2) EB > 80 %, the sample size of replication study can be determined.

Credible interval of the average power
The proposed interval estimation method can also be used to construct the credible interval of the average power. We resample the test statistics from the primary study N 1 times. In each run, a re-estimated variance of the effect sizes in the non-null SNPs σ 2 0 can be obtained. For a fixed σ 2 0 value, we first calculate the local true discovery rate of the primary study ltdr (1) with Eq. (17)
In simulation experiments, sample sizes are n 1 = 500 in primary study and replication study, respectively. The number of simulated SNPs is m = 1 × 10 4 . For each SNP, its minor allele frequency is uniformly distributed as U(0.05, 0.5). Their effect sizes are generated from the following distribution: For our hypothetical disease, its prevalence is 1 %. To test the marginal association between SNPs and the disease, log-odds ratio test is used. The significance levels are α 1 = 5 × 10 −5 and α 2 = 5 × 10 −3 in primary study and replication study, respectively. Figure 2 shows the histogram of the differences between observed effect sizes μ (1) and their underlying true values μ for identified true associations. We plotted separatively for the associations with positive effect (μ > 0) and negative effects (μ < 0). For the associations with pos-  (1) tend to be larger than the magnitude of their true effect. The observed effect sizes tend to be exaggerated itive effect, the mean value of the estimated effect sizes is larger than the mean value of true effect sizes. On the contrary, the values of the estimated effect sizes tend to be smaller than their true values for the associations with negative effect. In both of these two cases, the magnitude of the observed effect sizes tend to be exaggerated, which indicates that the winner's curse generally exists in the associations identified from primary study.
In order to check whether our EB based power estimator η (2) EB can overcome the winner's curse, the histogram of the differences between estimated values and true values is shown. As a comparison, we will show the corresponding histogram for power estimator by plugging in observed effect size first. A lot of methods have been proposed to overcome the winner's curse in terms of effect size estimation. CMLE and EB can be used directly in z-values of log-odds ratio test. The individual-level genotype data are also simulated so that bootstrap based bias reduction method BR2 [10] can also be used as a comparison (We modified the implementation code of BR2 so that log-odds ratio test can be used in the software). A direct thinking is to plug these adjusted estimators in the power calculation formula. The corresponding histograms for these three adjusted plug-in based estimators are shown as comparisons. In Fig. 3a, we use plug-in rule to estimate the replication study's power, where the observed effect size is plugged in. The estimated power is β (2) μ (1) . We plot the histogram of the difference between β (2) μ (1) and the true power values β (2) (μ) in the figure. The overestimated effect size makes the estimated replication study's power overestimated as well. Figure 3b plots the histogram of the difference between β (2) μ (1) CMLE and β (2) (μ). The winner's curse has disappeared, but there is a large downward bias in the estimated results. Equation (10) also introduces an Empirical Bayes estimator of the effect size, which reads μ (1) Figure 3c, d plot the histogram of β (2) μ (1) BR2 − β (2) (μ) and β (2) μ (1) EB − β (2) (μ), respectively. Large upward biases still exist in the histograms. In contrast, Fig. 3e plots the histogram of η (2) EB −β (2) (μ), where η (2) EB is our proposed EB based power estimator. The bias almost disappeared, indicating η (2) EB is better than other estimators of the replication study's power in terms of overcoming the winner's curse. The experiment has run 5 times, and the same conclusion holds in each run. The empirical biases of these five estimators can be seen in the Table 2. Table 3 shows the root mean square error (RMSE) of the five estimators β (2) μ (1) , β (2) and η (2) EB in the 5 runs. We can see that η (2) EB is better than other methods in terms of estimation accuracy.
To investigate the performance of the interval estimation, the coverage probability of 95 % credible intervals for all identified associated SNPs is shown in Table 4. The coverage probabilities by using pure Monte-carlo sampling are presented on the left side. In agreement with our analysis in the last section, the coverage probability is lower than the nominal value 95 %. The coverage probabilities by using modified method which combines Monte-carlo sampling and bootstrap are shown on the right side, which are closer to the nominal value. The credible interval given by the combined method is wellcalibrated. Figure 4 shows the average powers and their estimated results in the 5 runs. The true value of the average power is covered by 95 % credible interval in all runs. For comparison, we also show the estimated average power by using plug-in rule, i.e.β (2) μ (1) ,β (2) μ (1) CMLE , β (2) andβ (2) μ (1) EB , in each run. The estimated power for each association is also weighted by its local true discovery rate ltdr in each average power estimator. The figure shows that the EB based average power estimatorη (2) EB is much closer to the true valueβ (2) (μ). These results indicate that the estimated average power can be regarded as a proxy of the average power, which can be used to design replication study.
To check the performance of our method when the effect sizes of the associated SNPs are do not follow normal distribution, we also simulated data with the following distributed effect sizes: (20) and where t 5,0.02 is a scaled t distribution with degree of freedom 5 and scaling factor 0.2. The distribution of the associated SNPs' effect sizes follow the Gaussian mixture model in the second case. The average empirical biases and RMSE of all estimators in these two cases are shown in Tables 5 and 6, respectively. From the tables, we can see that our method is still better in terms of overcoming winner's curse and providing higher estimation accuracy.

WTCCC datasets
To give an application example of our proposed method, we will determine the sample size of replication study used for verifying the 6 human common diseases' findings from Wellcome Trust Case Control Consortium (WTCCC). The 6 diseases include coronary artery disease, Crohn's disease, hypertension, rheumatoid arthritis, type 1 diabetes and type 2 diabetes. Each disease has 2000 cases in the dataset of the primary study. There are 3000 shared controls among all datasets. The following quality control procedure is used in the primary study's datasets:  Histogram of β (2) (μ (1) )−β (2) (μ) β (2) (μ (1) )−β (2)  Histogram of β (2) (μ CMLE (1) )−β (2) Fig. 3 The histograms of the error in power estimation for all associated SNPs identified from primary study. a The replication study's power is estimated by plugging the observed effect size in power calculation formula, i.e. β (2) μ (1) . b The power is estimated by plugging the CMLE based corrected effect size μ (1) CMLE in power calculation formula, i.e. β (2) μ (1) CMLE . c The power is estimated by plugging in the BR2 estimator, i.e. β (2) μ (1) BR2 . d The power is estimated by plugging the EB based corrected effect size μ (1) EB in power calculation, i.e. β (2) μ (1) EB . e The power is estimated by EB based method, i.e. η (2) EB . The mean value of the estimation error is drawn with vertical dashed line. From the figure, it can be seen that η (2) EB has the smallest bias in power estimation. The biases for these 5 estimators are 0.144, -0.068, 0.045, 0.047 and 0.021, respectively Table 2 Empirical biases of power estimators of the replication study in the simulation experiments. The settings of the experiments can be seen in the main text β (2) μ (1) β (2) (2) (2) μ (1) , β (2) μ (1) CMLE , β (2) μ (1) BR2 and β (2) μ (1) EB are the plug-in based estimators by using observed effect size, CMLE, BR2 and EB in the effect size estimation. η (2) EB is proposed EB-based estimator. Bold face indicates the estimator achieving the smallest bias. In the experiments, η (2) EB behaves better than others in terms of bias reduction The significance levels used in primary study and replication study are α 1 = 5 × 10 −8 and α 2 = 5 × 10 −6 , respectively. The Control-to-Case ratio of the replication study is set to 1. The inferred hyperparameters π 0 and σ 2 0 can be seen in Table 7. With these parameters, the relationships between the estimated average power using EB based methodη (2) EB and sample size of replication study n (2) can be seen in Fig. 5. In conclusion, to achieve 80 % average power of the replication study, we will need 6885 individuals for coronary artery disease, 8092 individuals for Crohn's disease, 10,014 individuals for hypertension, 5291 individuals for rheumatoid arthritis, 4094 individuals for type 1 diabetes and 6988 individuals for type 2 diabetes. The detail about the sample sizes needed for different values of average power can be seen in Table 8.
As a comparison, we also plot relationships between β (2) μ (1) and n (2) in Fig. 5, whereβ (2) μ (1) is the estimated average power by plugging in observed effect Table 3 Root mean square error (RMSE) of power estimators of the replication study in the simulation experiments. The settings of the experiments can be seen in the main text β (2) μ (1) (2) μ (1) , β (2) μ (1) CMLE , β (2) μ (1) BR2 and β (2) μ (1) EB are the plug-in based estimators by using observed effect size, CMLE, BR2 and EB in the effect size estimation. η (2) EB is proposed EB-based estimator. Bold face indicates the estimator achieving the smallest RMSE. In the experiments, η (2) EB behaves better than others in terms of higher estimation accuracy sizes. For a given sample size n (2) , the estimated average power value using EB based method is much smaller than β (2) μ (1) . This is reasonable becauseβ (2) μ (1) is overestimated due to the winner's curse, which is alleviated in the EB based method. To achieve 80 % average power Fig. 4 The true value of average power is covered by the credible interval in simulation experiments. The EB based average power estimatorη (2) EB estimates more accurately than other plug-in based average power estimatorsβ (2) μ (1) ,β (2) μ (1) CMLE ,β (2) μ (1) BR2 and β (2) μ (1) EB . Average powerβ (2) (μ) is defined in the main text. A method combining Monte Carlo method and bootstrap is proposed to create the credible intervals of the replication study's average power. In the 5 runs of the simulation experiments, the underlying true values of the average power is covered by created intervals. Compared to other plug-in based average power estimators,η (2) EB is the closest toβ (2) (μ) Table 5 When effect sizes follow the distribution of Eq. (20), the average empirical bias and root mean square error (RMSE) of power estimators of the replication study in the simulation experiments Average β (2) μ (1) β (2) μ (1) CMLE β (2) μ (1) BR2 β (2) EB behaves better than others in terms of bias reduction and providing high estimation accuracy of the replication study, the sample size needed is 3023 for coronary artery disease, 3369 for Crohn's disease, 3788 for hypertension, 2748 for rheumatoid arthritis, 2706 for type 1 diabetes and 3095 for type 2 diabetes when usinḡ β (2) μ (1) as the estimator of average power. The sample sizes needed for other values of average power are listed in Table 9. These determined sample sizes are much smaller than the sample sizes determined by EB based method, indicating an underpowered study will be designed if we estimate power with observed effect sizes.
For coronary artery disease and type 2 diabetes, we obtained the publicly available summary statistics of the meta-analysis from two consortiums: CARDIoGRAM-plusC4D Consortium [22] and DIAGRAM Consortium [23], respectively. CARDIoGRAM GWAS is a metaanalysis of 22 GWAS studies of European descent involving 22,233 cases and 64,762 controls. The odds ratio calculated from high power CARDIoGRAM GWAS will be used as underlying true odds ratio to calculate the average power of the replication study for coronary artery disease in WTCCC. The average power obtained in this manner is denoted asβ (2) μ meta . Figure 5a plots the relationship betweenβ (2) μ meta and n (2) , which is the sample size needed in the replication study. The figure shows that our EB based power estimatorη (2) EB is very close to the power calculated using the results of CARDIo-GRAM GWAS. Also it can be shown thatβ (2) μ meta is in the credible interval we estimated. DIAGRAM GWAS is a meta-analysis consisting of 12,171 type 2 diabetes cases and 56,862 controls across 12 GWAS from European descent populations. Similar to CARDIoGRAM Table 6 When effect sizes follow the distribution of Eq. (21), the average empirical bias and root mean square error (RMSE) of power estimators of the replication study in the simulation experiments Average β (2) μ (1) β (2) μ (1) CMLE β (2) μ (1) BR2 β (2)  Bold face indicates the estimator achieving the smallest value in terms of bias or RMSE. η (2) EB behaves better than others in terms of bias reduction and providing high estimation accuracy in the experiments GWAS, the allele based odds ratio calculated from DIA-GRAM GWAS is used for calculating the average power of the replication study for type 2 diabetes in WTCCC. Figure 5f plots the relationship betweenβ (2) μ meta and n (2) . It can be shown that the result estimated by our EB based methodη (2) EB is close to the power calculated using the results of DIAGRAM GWAS.
If the values of the local true discovery rates ltdr (1) have nearly the same level for all identified associations in the primary study, the variance of the average power will be inversely proportional to the number of the associations. When the identified number is small in the primary study, the credible interval for the average power is rather wide. This can be illustrated in the study of hypertension, where there is only 1 association showing genome-wide significance. From Fig. 5c, we can see that the credible interval is rather wide. If we want to consider the credible interval for this situation, then the sample size can drastically increased.

Discussion
We propose to design replication study under the casecontrol setting where log-odds ratio test is used. The method can also be generalized to other tests within z-test scheme, such as regression slope test used for quantitative trait.
As described in [7], the winner's curse depends strongly on the power of primary study. For a high power primary study, most non-null SNPs will result in significant associations after random draws from the population. Hence, the bias will be small in this case. There are more and more high power studies conducted for common diseases by using pooling strategy or meta-analysis strategy, but the high power studies for rare diseases are limited. Hence, it is still helpful and necessary to propose a designing procedure for the replication study with the consideration of winner's curse.
With the development of the cost-effective sequencing technique, the targets of association studies extend from common variations to rare variants. A commonly used a b c d e f Fig. 5 The relationship between estimated average powerη (2) EB and the sample size of the replication study n (2) for 6 diseases of the WTCCC dataset: a coronary artery disease, b Crohn's disease, c hypertension, d rheumatoid arthritis, e type 1 diabetes, f type 2 diabetes. The Control-to-Case ratio of the replication study is set to 1. The significance levels used in the primary study and the replication study are α 1 = 5 × 10 −8 and α 2 = 5 × 10 −6 , respectively. As a comparison, the relationship betweenβ (2) μ (1) and n (2) are also shown in the figure. For a fixed n (2) ,β (2) μ (1) is much larger thanη (2) EB . In (a) and (f),β (2) μ meta is the average power estimator by plugging the log-odds ratio obtained from high power meta-analysis study: CARDIoGRAM GWAS [22] and DIAGRAM GWAS [23], respectively. It can be shown thatη (2) EB is close toβ (2) μ meta strategy to discover associations with rare variants is the collapsing method [24], in which several rare variants in a certain group are pooled together to enrich the signal. For each group, a "super variant" is constructed. If log-odds ratio test is adopted in testing the association between "super variant" and the disease, our method can be used directly for designing the replication study. Some limitations of our approach need to be mentioned.

1.
The assumption of our approach is that all SNPs' effect sizes are drawn independently from a two-component mixture distribution. Linkage disequilibrium widely exists in SNPs. Correlated genotype patterns can also introduce correlation between their effect sizes. The power estimation can be further improved by using correlation information in the prior set-up. 2. Our proposed method assumes the effect sizes of associated SNPs are normally distributed. This thin tail distribution may not be realistic. How to design of replication study with other heavy-tail prior needs to be discussed.