Theory

We now characterize the extra variance that is caused by relatedness within a given case-control study, and use this to compute the expected inflation factor δ. Starting from the definition of D, in Equation 2, we can write Var ^*[D] as

where i ≠ i′, j ≠ j′. We now need to determine how the value of this expression depends on cryptic relatedness.

Since G_i and H_j are Bernoulli trials, we have:

The following two terms in Equation 4 account for the possibility of departures from Hardy-Weinberg equilibrium in the sample. Assuming that these factors are independent of case-control status, we can write these as

where F measures the extent of the departure from Hardy-Weinberg equilibrium [4,19]. If, in fact, there is a different average level of inbreeding in cases than in controls [20], then we would replace F in Equation 7 and thereafter, with an average F across the cases and controls. (Notice that, unlike here, the inflation factor used by Devlin and Roeder was defined relative to the trend test, so that Hardy-Weinberg departures cancel out in their formulation.)

In our model, the controls are sampled randomly from the population. This means that the terms and are zero. This follows because, conditional on p, the fact that a random allele in the population is B, or b, provides no additional information about the genotype of another case or control in the sample. The assumption that controls are sampled randomly will usually be a good approximation, even if controls are specifically ascertained as not having the disease. As we will show below, the size of these covariance terms depends on the recurrence risk ratio for the phenotype, and the recurrence risk ratio for being unaffected is typically near one.

Next, since case alleles G_i are each similarly distributed, we can reduce Equation 4 by characterizing a single covariance between case alleles and then collecting the sum of all covariance terms that contain only case alleles. Given this, the Hardy-Weinberg equilibrium terms, and Equation 5, Equation 4 simplifies to:

where i ≠ i′. And now, finally, we need to evaluate under a model of cryptic relatedness. In order to do this, we first need to evaluate the probability that alleles in affected individuals share a common ancestor in generation t before the present. This will allow us to calculate the extra relatedness in cases due to the phenotype.

Recall that K_p is the population prevalence of the disease; K_t is the probability that a relative of an affected individual is also affected, given that the two individuals share a common ancestor t generations before the present; and that λ_t = K_t/K_p is the corresponding ratio of risks [15]. Next, let T _{(i,a)(i′,a′)} denote the coalescent time of allele copies a and a′ from individuals i and i′. In a slight abuse of notation, we will abbreviate T _{(i,a)(i′,a′)} as T_ii′. In what follows, individuals i and i′ are random (unphenotyped) draws from the population, except when specifically noted (e.g., Φ_i = aff indicates that i is affected). Then, using Bayes' rule, we can compute the coalescence rates for two chromosomes sampled from affected cases in the population as follows:

where P[T_ii′ = t] denotes the prior probability of coalescence in generation t, for random (unphenotyped) individuals. Next, using the assumption that affected and unaffected individuals coalesce with random chromosomes at the same rate (Equation 1), it follows that P[Φ_i = aff|T _ii′ = t] = K_p, and hence

Equation 9 produces a pleasingly simple result: the coalescence rate for chromosomes from affected individuals is increased by a factor that is closely related to the standard recurrence risk ratio.

Recurrence Risk for Relatives

The recurrence risk ratio is an important quantity in genetic epidemiology, and is widely measured [1]. For siblings, typical recurrence risk ratios for complex diseases range from around 2 to 50. For more distant relationships, the risk ratio declines approximately geometrically toward 1 as the number of meioses separating two relatives increases.

In our theoretical development, we will assume that disease inheritance is governed by a single additive gene [15], unlinked to the marker locus of interest. Other genetic models, including more complex models, behave similarly to this, except that the rate of decay of λ_t with increasing t may differ somewhat [15], leading to different coefficients in the cryptic relatedness term in Equation 16 below.

For the additive model, [15] obtained an expression for the recurrence risk ratio, λ_r, for any possible relationship, r, in terms of the recurrence risk ratio for full siblings, λ_s:

where Φ_r is the kinship coefficient between rth-degree relatives. For example, Φ_r = 1/4 between sibs, and decays by 1/2 for each increment to r. To connect λ_r to our model—which is written in terms of coalescent time t instead of r—we need to be more explicit about the mating patterns in the population model.

For example, under the standard Wright-Fisher model where individuals select their parents independently at random, most relatives are “half-relations”: half-siblings, half-first cousins, half-second cousins, etc. In that case, for t = 1, 2, 3, … , the corresponding kinship coefficients are Φ_r = 1/8,1/32,1/128, and so on. Then for example, for t = 2, λ_t − 1 = 4(λ_s − 1)/32. If instead, mating is purely monogamous, but partners are still chosen at random, then all relationships are “full”: full siblings, full cousins, etc. That is, for t = 1, 2, 3, … , the corresponding kinship coefficients are Φ_r = 1/4,1/16,1/64,….

In summary, λ_t may be much larger than 1 for the closest relatives, but it becomes approximately 1 if the common ancestor is more than just a few generations ago (> 10 or 15, say). This qualitative conclusion does not depend strongly on the assumed genetic model. Referring to Equation 9, this means that chromosomes from affected individuals have an excess probability of coalescing extremely rapidly (within the past few generations). If they do not, then they behave essentially like random chromosomes, for which coalescence takes place on timescales of thousands of generations in typical populations (Figure 1).

Open in a separate window

Figure 1

Coalescence Rates for Pairs of Random Chromosomes (Red) and for Pairs of Chromosomes from Affected Individuals (Green)

Notice that chromosomes from affected individuals have a small excess probability of coalescing very rapidly (i.e., in the most recent ten generations or so). Otherwise, their coalescence rates are essentially like those of random chromosomes. The region at the left-hand side of the graph between the red and green lines represents the excess probability of very recent coalescence among case chromosomes (denoted R in the text). This is what gives rise to the effect of cryptic relatedness. For larger t, the line for cases drops slightly below the line for random individuals, since both distributions integrate to 1. These plots assume an additive genetic model, with λ_s = 60, the “half”-relationships mating model, and a population size of 2,000. The line for cases was generated under the approximation that the excess relatedness is completely limited to the first n = 10 generations. In this case, the maximum coalescent probability for case chromosomes is 0.00275, when t = 1; R ≈ 0.00334. As expected, the mean coalescence time is ≈ 4,000 generations for both distributions. Alterations in n yield similar results (unpublished data).

The dynamics of this process are reminiscent of structured coalescent models with many demes [21–23]. In those models, two chromosomes from the same deme either coalesce with each other very quickly or escape into the population at large, and coalesce on a much longer time scale. These two phases have been described by John Wakeley as the “scattering phase” and the “collecting phase,” respectively [24]. An extreme example of this type of process (with selfing) was illustrated by Rousset [25].

Calculating the Inflation Factor

As described above, ancestral chromosomes of affected individuals coalesce at an increased rate during the most recent few generations (Figure 1), and otherwise behave essentially like random chromosomes. We now provide a heuristic derivation of the inflation factor δ; later we show that our expression closely approximates the results obtained in simulations. For simplicity, we consider the following approximation.

Let R be the excess probability of very recent coalescence for affected chromosomes relative to random chromosomes. That is,

where n might be taken as 10 or 15, say. Then write:

To evaluate , notice that there are two cases. With probability R, the two chromosomes coalesce very rapidly due to their shared phenotype. In that case, they share such a recent common ancestor that they are almost certainly identical by descent. In the second case, with probability 1 − R, the two chromosomes behave as random chromosomes, and their genotypes are independent Bernoulli draws from the population frequencies:

And finally, substituting Equations 11, 13, and 7 into Equation 3, we obtain

Equation 14 is worthy of discussion. When the simplest model of independence among sampled alleles holds, then δ = 1. The term containing F corresponds to Hardy-Weinberg departures, due to inbreeding for instance. The summation term corresponds to the effect of cryptic relatedness; the sum itself can be thought of as calculating the excess probability of identity by descent between chromosomes from affected individuals. Overall, the effect of cryptic relatedness increases linearly with the sample size m (for a given population size and λ_t).

Applications to Specific Models

In this section, we evaluate Equation 14 under a range of specific models, in order to determine when cryptic relatedness is likely to have a substantial impact on case-control studies. The models presented assume an additive genetic model, as described above. At first, we will assume that the population is of constant size N, so that the probability of coalescence in generation t, P[T_ii′ = t], is (1 − 1/2N)^t−1/(2N). After that, we turn to models with population growth. For simplicity, we set F = 0.

The Inflation Factor in Populations of Constant Size

Recall from Equation 10 that λ_t − 1 = 4Φ_r(λ_s − 1). Recall also, that when individuals select their parents independently at random, as in the standard Wright-Fisher model, that most relatives are “half-relations” (e.g., half-siblings, half-cousins, etc.), and then the kinship coefficients Φ_r are 1/8, 1/32, 1/128, … for t = 1, 2, 3, etc. Using δ _half to indicate this situation where individuals are related via “half-relationships,” it follows that

Noting that for small t (provided that N is not small), and that ∑ 2^−2t+1 converges quickly to 2/3; Equation 15 can be further approximated as

If instead, mating is purely monogamous, but partners are still chosen at random, then all relationships are “full”—e.g., full siblings, full cousins, etc., and the kinship coefficients are two-fold higher. The corresponding inflation factor, δ _full, is

indicating that the impact of cryptic relatedness is approximately doubled when there is fully monogamous pairing of parents, compared to when there is independent pairing of parents for each offspring.

Simulations

To check the accuracy of our analytical results, we generated population histories via Wright-Fisher simulation and estimated the inflation factor, δ, for a given disease and population genetic model, as described in the Materials and Methods section. Results are presented in Table 1, and compared to predicted results from Equation 16. The results show close agreement between the analytical prediction and the simulation results. In some cases, the analytical results slightly overestimate the inflation factor, probably due to the approximations used in relating Equation 9 to δ.

Table 1

Values of the Inflation Factor as a Function of Model Parameters, and a Comparison of the Simulated (_mean) and Analytical (δ_A) Predictions, for Populations of Constant Size

Open in a separate window

While the choice of an additive model for the phenotype (i.e., a heterozygote has exactly one-half the penetrance for the phenotype as a homozygote for the risk allele does) is mathematically convenient, alternative modes of inheritance (including multilocus models, or models with dominance components) are certainly likely in practice. Such models will have the impact of changing the rate of decay of λ_t, and hence the coefficient of the cryptic relatedness term in Equations 16 or 17. While we do not present a complete exploration of such models, we have performed a modest number of additional simulations under non-additive models. We have found that those results are qualitatively similar to the results presented above (unpublished data).

Intrinsic Constraints on δ

Table 1 shows the predicted impact of cryptic relatedness for a range of possible disease parameters. The magnitude of the inflation factor is fairly small for all parameter combinations shown, with a maximum value of 1.07. To make this more concrete, an inflation factor of 1.07 implies a quite modest excess of false positives: for instance, a fraction 1.5 × 10⁻³ of tests would be significant at the p = 10⁻³ level. As another example, consider a genetic model based loosely on a study of autism [26], where λ_s = 75, and K_p of 0.0004. Assuming the full-sibling model of relatedness, a sample size of 1,000, and a population size of 2.5 million (i.e., the number required to find that many cases), δ is just 1.02.

These examples notwithstanding, however, Equations 16 and 17 seem to suggest that δ can be made arbitrarily large simply by increasing the sample size m. But in fact, the space of sensible models is actually rather constrained. Since m cannot exceed K_p times the population size, there is a practical limit on m for a given λ_s and population size. Because of this constraint, it is difficult to construct biologically plausible parameter combinations that result in substantial inflation factors for randomly mating populations of constant size.

To be more specific, let K_s be the rate of disease in full siblings of an affected proband, i.e., K_s = λ_s K_p. Furthermore, let f be the fraction of all affected individuals in the population that are included in the sample. Then, noting that f = m/NK_p, Equation 17 can be rewritten as

Therefore, since f ≤ 1, for diseases where K_s is smaller than, say, about 1%, the inflation factor is negligible. The only way to get large values of δ is to have high values of K_s − K_p and nearly complete ascertainment of cases (high f). For instance, if K_s were 0.2 and λ_s were 4, then the inflation factor could be as large as 1.1, producing a small excess of false positives. But the latter calculation assumes complete sampling of affected individuals (f = 1), which would usually be difficult for a common disease.

In summary, in populations of constant size, the impact of cryptic relatedness is generally very small, unless (1) K_s is quite large—more than 0.2, say, and (2) f is near 1, meaning that there is nearly complete ascertainment of cases from the population. Hence, cryptic relatedness should not be a serious concern for most complex trait studies in stable populations, assuming random sampling of cases. As we will show in the next section, the situation is more serious for models with population growth.

The Inflation Factor with Changes in Population Size

We now consider a model that allows for changes in population size. Let N_t represent the population size at time t. Then, provided that the coalescent probability 1/2N_t is not especially large in any of the recent generations, and since λ_t − 1 decays as t increases, we can rewrite and simplify Equation 14 to

where again λ_t refers to the recurrence risk ratio for coalescence time t. Because (λ_t − 1) decays quickly toward zero, it is apparent that only changes in population size during the last few generations will impact δ. Moreover, for given values of m and λ_t, smaller population sizes in the past will produce higher inflation factors.

To check the accuracy of our results regarding demographic expansion, we modified the forward simulation procedure used above such that instead of a single N, we simulated exponential growth that began at time t _onset in the recent past starting at an initial population size N_A. For each subsequent generation t, the population size was determined by the equation N_{t +1} = N_t · e^α for a growth rate α such that the population size in the final generation is N_f. We performed at least 10,000 repetitions for each parameter combination, and the 95% standard error about the mean for each estimated δ was no greater than 0.01. In our analytic calculation, we assumed the “half” relationships model, as in Equation 15 and 16.

Results of the simulations, for a range of parameter values, are summarized in Table 2. Again, the theoretical prediction in Equation 19 is close to the simulated values. Under very recent growth models, can be substantial (as much as 2.5 for the extreme growth scenario shown). Under more realistic models of population growth, the effect of cryptic relatedness is smaller, but still non-trivial. Based on these results, it seems clear that the magnitude of growth is an important factor for determining δ. In populations that have grown rapidly from small size in the past few generations, cryptic relatedness may indeed lead to high inflation factors. It should be noted that many of the models presented have extreme growth; hence, the higher levels of cryptic relatedness shown here are likely to exceed anything seen in practice in human populations.

Table 2

Values of the Inflation Factor in Very Recently Expanded Populations

Open in a separate window

The qualitative difference between the equilibrium model and the population-growth model can be understood as follows. Consider two studies in which m affected individuals are sampled from each of two populations that have the same current size. If one population is of fixed size, while the other has grown rapidly from a smaller size, then the probability that two individuals are closely related is much higher in the growing population than in the equilibrium population. It follows from Equation 19 that this produces a higher inflation factor in the growing population than in the stable one.

Cryptic Relatedness with Biased Sampling

Thus far, we have considered models that assume “good” sampling design, in the sense that the sample of cases represents a random sample of the affected individuals in a population. We now consider the impact of sampling schemes that bias toward enrolling close relatives as cases in a study. For the previous models, we showed that with random ascertainment of cases, the inflation factor δ is maximized with complete ascertainment of cases from a population. The following models are instead motivated by the scenario in which a study enrolls only a small fraction of the affected individuals in a large population but, due to sampling biases, tends to recruit close relatives. Such situations might arise in practice if, for example, a patient at a clinic or in a study encouraged affected family members to visit the same clinic, or also to enroll in the study.

As an extreme, but simple example, consider first the situation in which the case sample consists of m(1 − σ) unrelated affected individuals, plus mσ/2 pairs of affected siblings (σ [0, 1]). The controls are all unrelated to anyone else. Assume furthermore that there is not inbreeding, so that F = 0 and the probability of recent identity-by-descent for chromosomes in siblings is 0.5. (For simplicity, we assume both in this and the next model that the sampling is from a sufficiently large population relative to m that we can approximately ignore the impact of cryptic relatedness apart from that induced by the biased sampling of siblings.) Then recall from Equation 14 that δ ≈ 1 + F + (m − 1)R where R is the (average) excess probability of recent coalescence, computed across all pairs of case chromosomes. In this model, a fraction σ/(m − 1) of the pairs of individuals are siblings. The probability that a randomly selected chromosome a in one sibling and a′ in the other sibling descend from the same parental chromosome is R = 1/4. Hence, for this model we obtain δ ≈ 1 + σ/4. At most, if the entire case sample is made up of sibling pairs, δ = 1.25. Any relatedness among the controls would further increase δ.

As a second simple example, suppose that a study recruits only a small fraction of affected individuals from a large population, but that recruits sometimes then encourage their siblings to enroll. Let the number of siblings of a recruited individual be Poisson with mean g, and let h be the probability that an affected sibling goes on to enroll in the study, independently for each affected sibling. Then the number of siblings of the initial recruit who enroll as patients in the study is Pois(ghK_s). After some algebra, it follows that the expected fraction of pairs of case individuals in the sample who are siblings is γ(γ + 2)/[(m − 1)(γ + 1)], where γ = ghK_s. Hence (again taking F = 0), we obtain

From these examples, it seems that biased sampling of cases can have a substantial effect on inflating the test statistics—though this is less dramatic perhaps than might have been expected. For example, suppose that index cases have an average of g = 2 siblings, that they refer affected siblings with probability h = 0.5, and that K_s = 0.4. Then the inflation factor δ ≈ 1.17.

Estimating coalescent probabilities in the Hutterites.

We estimated the coalescent probabilities for pairs of alleles in two individual Hutterites by the following. Starting from the affected individuals in the population, or from a matched random sample of individuals from the current population, we simulated the inheritance of a pair of randomly chosen chromosomes from different individuals, backward through time, from the present to the founders of the population. If the two chromosomes coalesced to a common ancestral chromosome within the pedigree, we counted the number of meioses back to that common ancestor, reporting the average number if the number of meioses was different on the two lineages. We repeated this procedure until we observed at least 500,000 coalescence events within the simulation. To estimate the mean inbreeding coefficient (F) in this sample, we used the same procedure as above except that we picked the two chromosomes from the same random individual, traced them backward in time, and determined how frequently those two chromosomes coalesced within the pedigree.

We thank Carole Ober for providing the marker, phenotypic, and genealogical data used for the Hutterite data analysis and for comments on the manuscript; Rebecca Anderson and Natasha Phillips for additional assistance in organizing and interpreting the phenotype data; and Catherine Bourgain, Graham Coop, William Wen, Sebastian Zöllner, and the anonymous reviewers for helpful comments or discussion. This work was supported in part by the National Institutes of Health (HG002772) and a Hitchings-Elion award from Burroughs Wellcome Fund to JKP; BFV received support from the above grant to JKP as well as NIH DK55889 to Nancy J. Cox and from a Genetics Regulation Training Grant NIH/NIGMS NRSA 5 T32 GM07197.