Interval estimation of genetic susceptibility for retrospective case-control studies.

BACKGROUND: This article describes classical and Bayesian interval estimation of genetic susceptibility based on random samples with pre-specified numbers of unrelated cases and controls.
RESULTS: Frequencies of genotypes in cases and controls can be estimated directly from retrospective case-control data. On the other hand, genetic susceptibility defined as the expected proportion of cases among individuals with a particular genotype depends on the population proportion of cases (prevalence). Given this design, prevalence is an external parameter and hence the susceptibility cannot be estimated based on only the observed data. Interval estimation of susceptibility that can incorporate uncertainty in prevalence values is explored from both classical and Bayesian perspective. Similarity between classical and Bayesian interval estimates in terms of frequentist coverage probabilities for this problem allows an appealing interpretation of classical intervals as bounds for genetic susceptibility. In addition, it is observed that both the asymptotic classical and Bayesian interval estimates have comparable average length. These interval estimates serve as a very good approximation to the "exact" (finite sample) Bayesian interval estimates. Extension from genotypic to allelic susceptibility intervals shows dependency on phenotype-induced deviations from Hardy-Weinberg equilibrium.
CONCLUSIONS: The suggested classical and Bayesian interval estimates appear to perform reasonably well. Generally, the use of exact Bayesian interval estimation method is recommended for genetic susceptibility, however the asymptotic classical and approximate Bayesian methods are adequate for sample sizes of at least 50 cases and controls.

Background

Association mapping of complex phenotypes in case-control samples involves analysis of tables of genotype/allele counts collected at a large number of genetic loci. Relating single-locus genotype and allele frequencies to the outcome is a basic analysis step even if complex interactions among loci are expected. Indeed, biologically realistic models that involve multiple interacting polymorphisms may induce considerable "marginal effects" associated with individual loci. Often the numbers of cases and controls are fixed in advance by the experimental design, and the multiple markers are typed. Then case/control proportions remain the same for all markers, whereas genotype and allele numbers in cases and controls are subject to the random sampling variation. The reverse is however of greater interest: what is the genetic susceptibility, or the probability that a random individual will have a particular outcome given that a particular genotype or an allele is observed at a locus? Unlike the odds ratio, this parameter is not invariant with respect to the prospective vs. retrospective sampling schemes. Therefore, the point estimate of this probability can only be obtained from genotype counts in cases and controls by assuming a particular value of the population prevalence. Another issue is the degree of uncertainty in the estimate, which is being investigated in this article via interval estimation. Classical, or "frequentist" confidence interval (CI) is a well established framework for interval estimation. Such an interval is a random quantity, and the probability statements are made about proportions of times a random CI covers the fixed population parameter. A more relevant question is often about the variability, or uncertainty associated with the estimate of the population value. In other words, we would like to be able to make statements about the lower and upper bounds for genetic susceptibility, thus interpreting the interval as fixed, and the susceptibility as random, where randomness is due to the limited amount of data about the parameter. Bayesian (e.g. "credible") intervals provide such interpretation, but sometimes are criticized for subjectivity associated with the choice of prior distributions. It is not generally the case that classical and Bayesian intervals should correspond to each other, however such correspondence is possible for certain combinations of the likelihood and prior distributions. For example, P. Altham [1] showed that when testing for the difference of two binomial proportions, frequentist Fisher's exact test can be viewed as a Bayesian test when assuming binomial likelihoods and Beta priors for the distributions of proportions. Thus, although no priors are explicitly assumed by Fisher's test, one can recover prior distributions that are indirectly implied by the test. Datta and Mukerjee [2] provide an extensive review of such matching probability problems and show that in most parametric cases suitable priors can be constructed that would match the Bayesian credible intervals of a given size to that of a frequentist confidence interval up to second order. As we show here, under a beta-binomial model there is a close connection between the two types of intervals for genetic susceptibility, which allows flexibility in the interpretation.

Results

Susceptibility intervals

Association studies are often "retrospective" in the sense that samples of cases and controls are determined by recruitment ("fixed"), and the genetic variants, e.g. genotypes AA and (not-AA) at a genetic marker, are random. Estimated genotype frequencies among cases or controls, e.g. are obtained directly from the following table as .

Assume AA is the risk genotype. If the samples were random with respect to both rows and columns, then the estimated susceptibility (penetrance, or the positive predictive value), φ = Pr( | AA), would simply be estimated by x/(x + y). Because the rows are fixed, we need the population prevalence, w. Then,

where w = Pr(), p = Pr(AA | ), and q = Pr(AA | ). If the rows are assumed to arise independently from binomial distributions, the maximum likelihood estimates of the parameters p and q are given by = x/n, and = y/m, respectively. The population prevalence, w, has to be estimated externally, rather than from the data provided by retrospective case-control samples. Let denote such an estimate of w. Hence, a point estimate of φ can be obtained as .

The expression given above involves the ratio of random variables, y/x, and hence obtaining the exact sampling variance is problematic. Approximate methods have been suggested for the ratio of binomial variables [3,4]. An approximate variance using the logit transformation allows us to obtain a simple expression as well as to separate the terms for p, q and w in the variance expression. To see this, define

Classical interval

The first-order Taylor series approximation gives the variance of as

where the first two terms refer to the variance of the relative risk of AA on the log scale. Pepe [5] presented extensive discussion on estimation of similar quantities in the context of biomedical research. The last term of (3) can be further approximated as

and requires the knowledge of which must come from an external source. When the range (r = w- w) of w is known, we may define , the range on the logit scale. Then can be approximated by (r*/6)2, since three standard deviations from the mean cover over 99% of the (centered) normal distribution, f (x, 0, σ), i.e. , where f(x, 0, σ) denotes the probability density function of a Normal distribution with mean 0 and standard deviation σ. Henceforth, unless otherwise mentioned we use . Browne [6] discussed this issue and similar approaches of relating the range to the standard deviation.

The CI for susceptibility is obtained by inverting the endpoints of the asymptotic normal interval for η,

e.g. the upper point is inverted as . As we know that , and , the estimated variance becomes

Approximate Bayesian interval

The variance given by (6) is infinite when the observed value of either x or y is zero. Before describing a method of dealing with this, we note that the usual asymptotic variance formula for the log odds ratios, based on the sum of reciprocals of the 2 × 2 table counts, {n11, n12, n21, n22}, has a similar deficiency, i.e. it results in an infinite variance when one of the observations is zero. To avoid this problem, it is common to add 1/2 to each cell. Haldane [7], Gart and Zweifel [8] justified this on the basis of minimizing the bias in the estimation of the logit variance (also see discussion in [9]). The variance becomes

We propose a similar modification of (6) and justify it by showing that it is an approximate Bayesian variance estimator. In passing, we note that (7) can also be obtained by using the same argument as in the derivation of the approximate posterior variance of that follows.

When the sampling distribution of genotypes is binomial, i.e. x|p ~ Bin(p, n) and y|q ~ Bin(q, m) and the prior distributions for p and q are independent Beta(γ 1, β1) and Beta(γ2, β2) respectively, the posterior distributions of p and q are independent Beta's with

p | x ~ Beta(x + γ1, n - x + β1)

q | y ~ Beta(y + γ2, m - y + β2)    (8)

The binomial likelihood arises as a consequence of independence between genotypes in a random sample. Although a Beta prior for allele frequencies is sometimes justified from purely mathematical convenience, it does provide a sufficiently flexible shape. More importantly, certain population-genetic models result in Beta, or more generally Dirichlet distributions for allele frequencies. For example, population frequencies at mutation-drift equilibrium follow this distribution [10]. A reasonable assumption is that a typical value from the susceptibility distribution is around the prevalence value, which implies similarity in the prior distributions of p and q, i.e. γ1 = γ2, β1 = β2.

Different parameters would generally assume prior deviation of the typical susceptibility value away from the population prevalence.

Posterior expectations and variances of p and q are respectively given by,

Using posterior variances for and , and substituting posterior expectations for parameter estimates and in (3), we obtain

Again, the approximation (4) is appropriate when is available.

When γ1, γ2 → 0 and β1, β2 → 0 (Haldane's prior, see [11]) we are essentially back to (6),

because n/(n + 1) and m/(m + 1) approach one as m and n increase. Therefore (6) is also justified from the Bayesian point of view.

When γ1 = γ2 = β1 = β2 = 1/2 (Jeffreys' prior, [11]),

Welch and Peers [12] established that in general Jeffreys' priors provide mathematical correspondence between frequentist and Bayesian intervals. In particular, these authors proved that under some regularity conditions (which are trivially satisfied for this problem) Jeffreys' prior is the unique prior for a parameter for which a Bayesian credible interval and a frequentist confidence interval match up to the second order.

The asymptotic normal credible interval based on can be obtained as:

where , and the endpoints, l and u, are inverted as and , respectively to produce an approximate (1 - α)% interval estimate of φ.

Exact Bayesian interval

It is also possible to obtain an "exact" Bayesian credible interval using Monte Carlo samples generated from the posterior distribution of φ. The sample is obtained by repeatedly drawing p and q from their posterior distributions given by (8) and calculating φ via (1) for each realization of p, q. This generates an empirical posterior distribution for φ and a (1 - α)% interval is given by (α/2), (1 - α/2) quantiles of this distribution. Values w can be drawn uniformly from the range of its possible values. A bell-shaped Beta distribution can also be assumed with parameters based on the reported range, w- w, as described in Methods, section 3. An algorithm is as follows:

1. Generate p(~ Beta(x + γ1, n - x + β1), i = 1,.., B.

2. Generate q(~ Beta(y + γ2, m - y + β2), i = 1,.., B.

3. Generate w(~ Uniform (w, w), i = 1,.., B.

4. Compute , i = 1,..., B.

Then the α/2 and 1 - α/2 quantiles of the generated empirical posterior distribution of φ provide the (1 - α)% credible interval for the susceptibility. This algorithm is easily implemented in R, a language and environment for statistical computing and graphics available as Free Software at . For example, setting γ1 = γ2 = β1 = β2 = 0.5, the entire code using 100,000 samples, 95% credible interval, and w distributed uniformly between 0.04 and 0.06 is

p <- rbeta(100000, x + 0.5, n - x + 0.5)

q <- rbeta(100000, y + 0.5, m - y + 0.5)

w <- runif(100000, 0.04, 0.06)

phi <- w*p / (w*p + (1 - w)*q)

phiL <- quantile(phi, prob = 0.025)

phiU <- quantile(phi, prob = 0.975)

cat(c("Interval estimate", phiL, phiU), fill = T)

Common beta priors on p and q (γ1 = γ2, β1 = β2) may be interpreted as vague with respect to the prior distribution for φ. On the logit scale, this distribution is bell-shaped symmetric, and centered around ln [w/(1 - w)]. Its variance decreases as the common (γ, β) parameters for p and q increase. For large values of γ, β this variance can be obtained by the Taylor series approximation, 2β/[γ (γ + β + 1)]. A mixture of beta distributions can take an arbitrary shape with the great advantage that the resulting distribution is bounded between the lower and upper values – which is the case for p, q, and φ parameters. In Methods (section 1) we outline the sampling scheme for the mixture.

Interval estimation of allele or haplotype susceptibilities

Methods described so far could be applied without modifications to the estimation of allele or haplotype susceptibilities when the probability of obtaining a sample of alleles follows the multinomial likelihood. However, such usage would require the assumption of independence between alleles. This can be justified under the assumptions of Hardy-Weinberg equilibrium (HWE) in the population, as well as the multiplicative effects of haplotype susceptibilities [13]. When these assumptions are reasonable, the intervals are obtained in the same way as described above using counts of alleles/haplotypes instead of genotypes. Otherwise the variance of allele frequencies, is no longer binomial (x, y are now the counts of alleles A in cases and controls respectively, and n, m are the total numbers of alleles). It is instead given by

where D, Dare the Hardy-Weinberg disequilibrium (HWD) coefficients [14], defined as the deviation of the observed frequencies of AA from the expected under random union of alleles. For example, , and is estimated by using the observed frequencies, . These variances are most easily derived by re-coding the genotypes (AA, , ) as G = (-1, 0, 1), which keeps track of the number of alleles A minus one in a corresponding genotype, and taking expectations in V(G) = E(G2) - (E(G))2 (see Appendix 1 in [15], or two alternative derivations in [14]).

These coefficients (quantifying HWD among cases or controls) may be non-zero even if the population is in HWE [13]. Then the variance (6) becomes

with two new additional terms to account for HWD. Consider the effect of . If is negative, the resulting variance becomes smaller. If it is positive, the maximum possible value of [14] is

The maximum values for the terms with HWD in (15) are

Therefore, when or reach their maximum possible values, the inflation of the corresponding part of the variance is twice the value assuming the equilibrium. To set up an exact Bayesian interval for allelic/haplotypic susceptibility we may start with distributions for genotypes/diplotypes instead of that for alleles/haplotypes. Random samples of genotypes for cases, {n1}, and for controls, {n2}, follow the multinomial distribution. A conjugate prior distribution is Dirichlet(γ), so the posterior genotype frequencies Pare sampled from Dirichlet(γ+ n) distribution. Posterior allele/haplotype frequencies are given as and the method proceeds as above. Jeffrey's prior is obtained when all γare 1/2.

Coverage probabilities

The intervals using (6, 10) and the exact Bayesian interval have good frequentist coverage probabilities. To illustrate this, we simulated population values of p and q as p ~ Beta(1, 1) and q ~ Beta(1, 1) for each of 10,000 simulation runs. We used γ = β = 1/2 for the interval calculations and assumed that w = 0.5 or w = 0.04 is known and fixed. To compare intervals based on (6) and (12) we only considered samples with x, y > 0, obtained as x ~ Bin(p, n), y ~ Bin(q, m). In addition, we considered binomial samples with at least one of x, y equal to zero to check coverage properties of (12). We used m = n = 10,50,100,500 and 90% coverage intervals. Frequentist intervals cannot be used with zero counts of n11 or n21, which may be problematic only for sample sizes of 10 (17% of the samples). In this small percentage of cases, we resort to adding 1/2 to the counts as in (10). Agresti [9] considered a similar approach for calculation of CIs for odds ratios, that is to use Gart's formula (7) only when zero counts are encountered. Results for coverage probabilities are shown in Table 1. Note that the coverage of all three methods is around the nominal 90% for all sample sizes and both values of population prevalence. Results in Table 1 are averages across the distribution of p and q, however similar coverage results are obtained when the population values of p, q are set to specific fixed values. We considered the following fixed sets of parameters:

Table 1

Coverage probabilities for nominal 90% intervals based on 10,000 simulations.

Sample size	10	50	100	500
Prevalence
	Asymptotic Frequentist
0.04	0.921	0.910	0.908	0.898
0.50	0.918	0.906	0.906	0.898
	Asymptotic Bayesian
0.04	0.902	0.903	0.904	0.897
0.50	0.896	0.898	0.901	0.898
	Exact Bayesian
0.04	0.909	0.902	0.902	0.896
0.50	0.907	0.895	0.901	0.895

1. p = 0.96, q = 0.06, w = 0.04. This set results (equation 1) in the population susceptibility value of φ = 0.4 (ten-fold increase from the population prevalence value).

2. p = 0.9, q = 0.15, w = 0.04. This set results in the population susceptibility value of φ = 0.2.

3. p = 0.8, q = 0.2, w = 0.5. This set results in the population susceptibility value of φ = 0.8.

The coverage probabilities for all three methods were found to be sufficiently close to the nominal values of 90% and 95% for sample sizes of 50 or larger (Table 2).

Table 2

Coverage probabilities for nominal (90%/95%) intervals based on 10,000 simulations, based on three population settings of (p, q, w). Setting 1: p = 0.96, q = 0.06, w = 0.04. Setting 2: p = 0.9, q = 0.15, w = 0.04. Setting 3: p = 0.8, q = 0.2, w = 0.5.

Sample size	10	50	100	500
(p, q, w) settings
	Asymptotic Frequentist
1	0.888/0.951	0.919/0.962	0.916/0.960	0.896/0.948
2	0.942/0.883	0.906/0.962	0.909/0.954	0.896/0.953
3	0.909/0.954	0.908/0.954	0.900/0.947	0.899/0.951
	Asymptotic Bayesian
1	0.888/0.883	0.918/0.920	0.907/0.951	0.901/0.948
2	0.866/0.920	0.900/0.950	0.904/0.946	0.894/0.951
3	0.879/0.907	0.901/0.950	0.895/0.944	0.899/0.948
	Exact Bayesian
1	0.911/0.979	0.890/0.922	0.894/0.954	0.897/0.946
2	0.945/0.967	0.900/0.950	0.908/0.950	0.896/0.951
3	0.865/0.954	0.900/0.948	0.896/0.946	0.895/0.950

Table 3 shows the mean length and the standard deviation of the interval length across 10,000 simulations used to produce Table 1. This standard deviation reflects estimated variability of the interval length and decreases with sample size, but not with the number of simulations. The corresponding standard error could be obtained by dividing the standard deviation by the square root of the number of simulations. Table 4 shows the same results among the intervals that include the population value of w. Numbers for each parameter combination that contrast three intervals in a given table are obtained using the same data sets. For example, once a sample of size n is obtained for a specific value of w, all three intervals are calculated using this sample. Therefore, comparisons between methods are quite precise given the number of simulations. Again, the mean and standard deviation are similar among three methods. Nevertheless, the classical CI is never the shortest of the three intervals for the parameter values that we have considered and typically it has the highest length variability. Using data summarized in Table 3, the exact interval was found to have the smallest length 5 times out of 8 considered parameter combinations and has the lowest variability 5 times out of 8. These numbers are respectively 5 out of 8 and 4 out of 8 for data used to construct Table 4. There is some dependency between results from Tables 3 and 4 due to the same combination of parameters.

Table 3

Average length (standard deviation) of the intervals based on 10,000 simulations.

Sample size	10	50	100	500
Prevalence
	Asymptotic Frequentist
0.04	0.161(0.220)	0.081(0.151)	0.060(0.121)	0.026(0.057)
0.50	0.361(0.156)	0.168(0.090)	0.117(0.062)	0.052(0.030)
	Asymptotic Bayesian
0.04	0.143(0.207)	0.075(0.139)	0.057(0.114)	0.025(0.055)
0.50	0.341(0.146)	0.166(0.086)	0.116(0.060)	0.052(0.030)
	Exact Bayesian
0.04	0.180(0.264)	0.079(0.146)	0.058(0.113)	0.025(0.054)
0.50	0.329(0.152)	0.161(0.085)	0.114(0.059)	0.052(0.029)

Table 4

Average length (standard deviation) of those intervals that contained the true parameter value.

Sample size	10	50	100	500
Prevalence
	Asymptotic Frequentist
0.04	0.157(0.215)	0.078(0.147)	0.059(0.117)	0.026(0.056)
0.50	0.364(0.155)	0.168(0.088)	0.119(0.067)	0.052(0.030)
	Asymptotic Bayesian
0.04	0.139(0.201)	0.072(0.134)	0.056(0.110)	0.026(0.055)
0.50	0.342(0.144)	0.165(0.084)	0.118(0.065)	0.052(0.030)
	Exact Bayesian
0.04	0.176(0.259)	0.077(0.142)	0.057(0.111)	0.025(0.053)
0.50	0.333(0.150)	0.161(0.084)	0.116(0.065)	0.052(0.030)

We do not present numerical results based on mis-specification of w or the uncertainty in it, because the behavior of intervals is clear. When w is mis-specified, the intervals actively worsen as the sample size increases (in the sense of the probability of including the population parameter), since they shrink around the wrong value. On the other hand, the uncertainty in w widens the interval length. As the sample size becomes large, the minimum length simply reflects the variability in w. We conducted simulation experiments to illustrate this point, using the population value of w = 0.045 assuming the range 0.03 to 0.06 is reported. For the exact Bayesian interval, the distribution of w was modelled as Beta(ψ1 = 77.31, ψ2 = 1640.69) with parameters estimated by (29). These values are taken to reflect the reported prevalence of the hypersensitivity reaction to an antiretroviral drug abacavir, for which genetic predisposition has been described [16]. Table 5 reports results for the coverage probabilities and the interval length of 10% intervals, assuming p = 0.8 and q = 0.2, which corresponds to the population susceptibility value of 0.159. The results show that the average length is somewhat increased when the calculation are based on the range of w, and the interval coverage is above the nominal value of 90%, because of the increased interval width.

Table 5

Effect of uncertainty in the prevalence, w: the population value is w = 0.045; the assumed range is either 0.03–0.06 ( > 0, first three columns), or zero ( = 0, last three columns). AF, AB, EB refer to the asymptotic frequentist, approximate Bayesian, and exact Bayesian 10% intervals.

	> 0			= 0
Interval	AF	AB	EB	AF	AB	EB
n = m = 10
Av. Length	0.375	0.308	0.394	0.371	0.303	0.391
Coverage	0.941	0.889	0.868	0.913	0.889	0.868
n = m = 50
Av. Length	0.154	0.144	0.155	0.143	0.134	0.145
Coverage	0.929	0.924	0.920	0.902	0.898	0.895
n = m = 100
Av. Length	0.110	0.107	0.110	0.096	0.093	0.097
Coverage	0.945	0.946	0.939	0.900	0.896	0.900
n = m = 500
Av. Length	0.067	0.067	0.066	0.041	0.041	0.041
Coverage	0.993	0.993	0.992	0.897	0.899	0.900

Pharmacogenetic application

Hypersensitivity reaction to the antiretroviral drug abacavir affects approximately 4.5% of patients. The hypersensitivity symptoms are varied with fever and rash among the most common. A small number of fatalities have also been reported. The exact mechanism of hypersensitivity is not established, but the pattern of symptoms suggests that it is an immunological reaction, triggered by specific genetic polymorphisms. Hetherington et al. [16] studied the association of HLA-B57 haplotype with hypersensitivity. One or two copies of HLA-B57 haplotype were present in x = 39 of n = 84 individuals with hypersensitivity (cases) and in y = 4 out of m = 113 controls. The prevalence of the hypersensitivity reaction is estimated to be between 0.03 and 0.06 [16]. The estimate of φ using (1) and the midrange w (that is equal to 0.045) is 0.382. Assuming w = 0.045 is fixed (i.e. V(w) = 0), the intervals are as follows. The frequentist 95% interval is 0.187 - 0.624, the asymptotic Bayesian interval is 0.180 - 0.584, and the exact Bayesian interval is 0.203 - 0.648. Taking the range of w into account, the frequentist and Bayesian asymptotic intervals become wider, 0.183 - 0.631 and 0.175 - 0.591, correspondingly, and the exact Bayesian interval assuming the uniformity of w on the range (0.03 - 0.06) is 0.186 - 0.658. This interval is estimated from repeatedly (i = 1,..., 100000) generating p(~ Beta(39 + 1/2, 84 - 39 + 1/2), q(~ Beta(4 + 1/2, 113 - 4 + 1/2), w(~ Uniform(0.03,0.06), and calculating . Then the credible interval was obtained from this (empirical) distribution of φ by determining its 2.5% and 97.5% quantiles. It is not necessary to assume the uniformity of w. Indeed, it may be reasonable to assume that the midrange value is more plausible than the endpoints and a realistic distribution can be described as bell-shaped. Taking w ~ Beta(57.8572,1227.86) generates a bell-shaped distribution with the mean at 0.045 and 99% of this distribution is contained between 0.03 and 0.06 (Methods, section 3). The resulting interval is somewhat smaller, 0.195 - 0.653. Results for this example are summarized in Table 6.

Table 6

Summary of susceptibility intervals for the pharmacogenetic example. AF, AB, EB refer to the 95% asymptotic frequentist, approximate Bayesian, and exact Bayesian intervals.

	AF	AB	EB
= 0	0.187–0.624	0.180–0.584	0.203–0.648
> 0	0.183–0.631	0.175–0.591	0.195–0.653

Discussion

Asymptotic frequentist and Bayesian intervals for genetic susceptibility as well as the exact Bayesian interval are quite similar in properties. Given that the sampling from the posterior distribution can be done directly, there is no particular reason to resort to approximations and the exact Bayesian interval can be recommended. Nevertheless, approximations are useful in that they reveal connections between the confidence and credible intervals as well as allow for sample size and power calculations as we outline in Methods, section 2. The mean length and the interval variability appear somewhat smaller for the Bayesian intervals. The frequentist coverage of Bayesian intervals is satisfactory. One may even argue that the operational usage of CIs is more often Bayesian than it is not, because the practical issue is ultimately about the confidence on the plausible range of parameter values. Algebraic relations between classical and Bayesian intervals as well as similarities in numerical properties allow Bayesian interpretation for the CIs. On the other hand, Bayesian intervals for this problem can be described in the classical sense as random, covering the fixed population parameter α% of the time in repeated samples.

Caution should be taken when the inference is about allele or haplotype susceptibilities. When assumption of population HWE does not hold, or when there are substantial deviations from multiplicativity, the variance of susceptibility contains additional terms that account for population or model-induced deviations from Hardy-Weinberg equilibrium.

Conclusions

We found that the classical interval for genetic susceptibility can be considered as an approximate Bayesian interval under the beta-binomial model. This algebraic similarity between the proposed classical and approximate Bayesian intervals allows an appealing Bayesian interpretation of the usual confidence limits. Simulation studies also confirm similarities in coverage probabilities and interval lengths among asymptotic classical and approximate and exact Bayesian intervals.

Methods

Posterior sampling using the mixture of beta distributions

We derive explicitly the posterior distribution of θ under a Binomial sampling using a mixture of Beta priors. As Beta is a conjugate prior for the binomial likelihood, we show that mixture of Beta's is also conjugate. To see this, consider the problem

where p(x| θ) denotes the sampling density of X = x given θ. Now consider a prior for θ given by a mixture of k beta distributions Beta(γ, β) with weights v, such that and v≥ 0 for all j. More explicitly, the prior density π(θ), of θ is given by,

where denotes the density of a Beta (γ, β) distribution and denotes the Beta function. Using Bayes' rule it follows that the posterior distribution is given by,

where = γ+ x, = β+ n - x and

Thus it follows that if X | θ ~ Bin(θ, n) and θ ~ vBeta(γ, β) then θ|X = x ~ . It may be noted that any continuous density on [0,1] can be approximated by an appropriate mixture of Beta densities (see [17] for more general results on approximating any prior using natural conjugate priors). Besides such theoretical properties, sampling from the mixture of Beta's is also straightforward. For instance, in order to generate θ ~ , we can apply the following two-stage sampling:

1. Sample J = j ∈ {1,..., k} with probability , i.e., Pr [J = j] = .

2. Conditional on the sampled value of J = j, sample

This generates a sample from the mixture . Thus, our sampling scheme described in "Exact Bayesian Interval" section can be easily modified to sampling from mixture of Beta's instead.

Power and sample size

Under the hypothesis of equality of genotype frequencies in cases and controls, H0 : Pr(AA | ) = Pr(AA | ), the expected susceptibility in (1) is equal to the population prevalence, w. Correspondingly, the logit transformation η has the mean equal to , and

The two-sided alternative hypothesis is H: Pr(AA | ) ≠ Pr(AA | ). The values of Pr(AA | ), Pr(AA | ) and the ratio under the alternative are denoted as p, q, and Δ, respectively. Once these values are specified, we can calculate the power of detecting the difference of the susceptibility from the population prevalence as well as the sample size required to achieve the required power under an acceptable type I error rate. We will illustrate the sample size calculation assuming the equal number of cases and controls (m = n). Power calculation given the fixed sample size is obtained similarly. We have

where

is the non-centrality parameter corresponding to 1 - β power under a two-sided level-α test. Therefore, the necessary sample size is

Estimating parameters of the Beta distribution for the prevalence from the range of reported values

Suppose the range of the prevalence values, w- w, is reported. We would like to model the uncertainty with a bell-shaped Beta(ψ1 > 1, ψ2 > 1) distribution, such that a certain percentage of the distribution (e.g. 99%) is between wand wwith the expected value given by w+ (w- w)/2. The parameters of this distribution can be found as the solution for ψ1 of

where denotes the density of a Beta(ψ1, ψ2) distribution and is the Beta function. The second parameter for b(·) in (28) is given by the constraint that the expectation for w is ψ1/(ψ1 + ψ2).

As a sample calculation, consider an antiretroviral drug abacavir for which the hyper-sensitivity reaction range is reported as w= 0.03 to w= 0.06 [16], so that we can take = w+ (w- w)/2 = 0.045. If we assume that these bounds cover 99% of the prevalence distribution that results from sampling errors in estimation, possibly confounded with other factors such as genuine population heterogeneity of samples, then the resulting distribution reflecting the uncertainty in w is Beta(57.8572, 1227.86), using (28).

Much more simply calculated parameter estimates can be obtained by the normal approximation to the beta distribution. Let = w+ (w- w)/2 and approximate the variance by = [(w- w)/6]2. By the method of moments, estimates of the parameters (ψ1, ψ2) are

For the abacavir example, these estimates are found as = 77.31 and = 1640.69. Using these estimates, 99.7% of the resulting Beta(77.31, 1640.69) distribution is contained within the range 0.03 - 0.06. Thus, the normal approximation yields a somewhat more condensed distribution.

Authors' contributions

All authors made substantial contributions to this paper, including conceiving of the ideas, discussion and writing. All authors read and approved the final manuscript.