Abstract
Relatedness between individuals is central to many studies in genetics and population biology. A variety of estimators have been developed to enable molecular marker data to quantify relatedness. Despite this, no effort has been given to characterize the traditional maximumlikelihood estimator in relation to the remainder. This article quantifies its statistical performance under a range of biologically relevant sampling conditions. Under the same range of conditions, the statistical performance of five other commonly used estimators of relatedness is quantified. Comparison among these estimators indicates that the traditional maximumlikelihood estimator exhibits a lower standard error under essentially all conditions. Only for very large amounts of genetic information do most of the other estimators approach the likelihood estimator. However, the likelihood estimator is more biased than any of the others, especially when the amount of genetic information is low or the actual relationship being estimated is near the boundary of the parameter space. Even under these conditions, the amount of bias can be greatly reduced, potentially to biologically irrelevant levels, with suitable genetic sampling. Additionally, the likelihood estimator generally exhibits the lowest root meansquare error, an indication that the bias in fact is quite small. Alternative estimators restricted to yield only biologically interpretable estimates exhibit lower standard errors and greater bias than do unrestricted ones, but generally do not improve over the maximumlikelihood estimator and in some cases exhibit even greater bias. Although some nonlikelihood estimators exhibit better performance with respect to specific metrics under some conditions, none approach the high level of performance exhibited by the likelihood estimator across all conditions and all metrics of performance.
AN understanding of the relatedness between individuals plays an important role in many areas of population biology and genetics. For example, it is central to quantitative genetics and plays a crucial role in estimating heritability and additive genetic variances and covariances (Falconer 1981; Lynch and Walsh 1998). Likewise, it may be useful in studies of isolationbydistance or population structure. Consequently, a number of different means of quantifying relatedness have been developed. Most inclusive of these are the sets of identitybydescent modes described by Jacquard (1974). However, for large noninbred populations these reduce to a pair of quantities: the probability that two individuals share two alleles identicalbydescent and the probability that they share one allele identicalbydescent. More commonly, the coefficient of coancestry θ (Jacquard 1974) or the coefficient of relatedness r = 2θ are used to quantify the degree of relatedness between two individuals.
Estimates of θ or r may be derived in a variety of ways. Traditionally, they are calculated from a known pedigree (Crow and Kimura 1970). Increasingly, however, molecular marker data have been used to estimate relatedness. Consequently, a number of estimators have been developed for this purpose (Thompson 1975; Queller and Goodnight 1989; Liet al. 1993; Ritland 1996a; Lynch and Ritland 1999; Wang 2002). These have been developed in a variety of different ways. The Queller and Goodnight (1989) estimator was motivated to ensure that Hamilton’s rule (Hamilton 1964a,b) applied under general circumstances given an estimate of relatedness r; in contrast, the Ritland (1996a), Lynch and Ritland (1999), and Wang (2002) estimators were based on different methodofmoments approaches to the relationship between relatedness and genotypic similarity.
Both Ritland (1996a) and Lynch and Ritland (1999) mention maximumlikelihood estimators of relatedness coefficients; however, they dismiss their utility on the basis of simulations that indicate that many (e.g., 70 or more) loci may be required. However, both approaches deviate from the traditional approach involving a likelihood function defined for the set of three parameters sufficient for describing relatedness in a noninbred population (Thompson 1975). While the likelihood function used by Ritland (1996a, p. 180) is a special case applicable to a single twogene relatedness parameter (e.g., θ), that proposed by Lynch and Ritland (1999, Equation 12) cannot be derived from the traditional one (see appendix b). Further, both admitted solutions outside the biologically meaningful parameter space (Thompson 1975) and thus yield estimates that cannot be interpreted as probabilities of identitybydescent. Thus, despite extensive efforts to characterize the statistical behavior of other estimators of relatedness, it remains unclear how the traditionallikelihood one compares.
The statistical behavior of relatedness estimators is critical to their utility in practice. Because of their complexity, simulations have generally been relied upon to characterize the sampling error. These simulations have taken two approaches. The first constructed data sets from relatively simple conditions, assuming identical allele frequency distributions across loci (Ritland 1996a; Lynch and Ritland 1999). A second approach (Van de Casteeleet al. 2001) was motivated by actual microsatellite data. In this case, data sets were constructed from more complex situations involving variation among loci in allele frequency distributions and variation in population structure. A similar approach involving actual data and a known pedigree was used to test the utility of markerbased estimators (Thomaset al. 2002). These two basic approaches have yielded somewhat contradictory pictures of the range of estimators available. It is clear, however, that not all estimators perform equally well under all conditions. Further, the same estimator may not perform best under all conditions. Which estimator performs best may depend not only on the nature of the biological conditions, but also on the criterion used to measure performance.
Conspicuously lacking from the array of relatedness estimators under test is the traditional general maximumlikelihood estimator (Thompson 1975). Often maximumlikelihood estimators exhibit many desirable features (Kendallet al. 1979), including having lower standard error, being asymptotically unbiased, being adaptable to a wide variety of sampling conditions, and naturally accounting for differences among different subsets of the sample, for example, different allelefrequency distributions among loci. The primary negative feature of maximumlikelihood estimators is that they are often biased for finite sample sizes. However, in many cases that bias is small enough to be biologically irrelevant or can be dramatically reduced with reasonable sampling.
Given the many desirable features of maximumlikelihood estimators, it would be natural to develop one for relatedness of individuals. Further, it would be useful to determine whether a maximumlikelihood estimator of relatedness approaches its asymptotic properties rapidly enough to be useful in practice or to compete with nonlikelihood estimators. In this study, we investigate one such estimator (Thompson 1975) on the basis of the genetic information available for two individuals assayed at many different loci. In addition to the genotypic information, the estimator relies on the allele frequency distributions at each locus sampled. For this analysis we assume that the allelefrequency distributions are known without error to focus on the behavior of the relatedness estimator itself. In practice, realistic samples will often involve enough individuals that errors in the allelefrequency distribution will be quite small. In some cases, for example, microsatellites segregating many alleles, errors in the frequencies may be significant. The importance of the additional sampling variance introduced is beyond the scope of this article.
Note that the approach taken here, to estimate relationship as a continuous parameter, is distinct from a related one, which infers the degree of relatedness from among a set of discrete possibilities. Both approaches are discussed by Thompson (1975). The latter approach is further developed by Thompson (1986) and has been used extensively in human genetics recently to classify groups of individuals (usually pairs) into distinct relationship classes (Boehnke and Cox 1997; Painter 1997; Broman and Weber 1998; Siebertset al. 2002).
The primary goal of this study is to assess the performance of the likelihood estimator, in comparison with some of those already developed, under a range of biological sampling conditions. In particular, the performance is quantified by two measures of the distribution of estimates obtained for each estimator (the standard error and the bias) and by the overall deviation of those estimates from the parametric value, quantified by the root meansquare error. From this information it is possible to determine how aspects of the sampling conditions, e.g., number of loci or segregating alleles, or aspects of the relationship being estimated influence the ability of one estimator or another to perform well. Thus, some guidance for experimental design can be developed.
The initial focus of this study is on a relatively simple set of conditions, although one meant to mimic a variety of natural situations. In this sense, it is more closely connected to the evaluation used by Lynch and Ritland (1999) than to that used by Van de Casteele et al. (2001). This approach has been chosen to provide a clearer indication of how each different factor, e.g., number of loci sampled, number of segregating alleles, and allele frequencies, influences the performance of each estimator.
STATISTICAL MODELS
For population samples in which the ancestry of individual alleles is unknown, the most general means of describing the relatedness of one individual to another is in terms of the nine identity modes described by Jacquard (1974) (Figure 1). The degree of relatedness is quantified by a set of coefficients Δ= (Δ_{1}, Δ_{2},... Δ_{9}), each of which represents the probability of the four alleles at a single locus in two diploid individuals sharing the corresponding particular pattern of identitybydescent. In a large, noninbred population, only Δ_{7}, Δ_{8}, and Δ_{9} are nonzero; consequently, in such a population any pattern of relationship between two individuals can be described by that set of three coefficients. The most commonly used summary of the degree of relationship is the coefficient of coancestry (Jacquard 1974; Lynch and Walsh 1998),
Likelihood models: Likelihood estimators are based on a probability model of the sampled data. In this case, the unit of sampling is a pair of individuals, each one of which has been assayed genetically at L loci. The estimator described here is based on the assumption of independently segregating marker loci. The likelihood for the overall sample, therefore, is simply the product of the likelihoods across the loci.
The basic probability model of the sampled alleles at a single diploid locus is well known (Thompson 1975). Usually it is given only for the case of large, noninbred populations where only three modes of relatedness (S_{7}, S_{8}, and S_{9}) are possible. However, the structure of the model is much clearer in its general form and it is applicable to the full range of population structures. As a basis of further generalization it thus warrants explicitly outlining the complete onelocus likelihood model for the nine relatedness coefficients.
There are nine distinct patterns of identityinstate for the four alleles sampled at a single locus in two individuals. Table 1 lists these in the second column. As an example, a pair of individuals each homozygous for allele A_{1} represent identityinstate mode
Given that a pair of individuals is known to be related according to identitybydescent mode S_{j}, the probability of each identityinstate pattern S_{i} is dependent on the allele frequencies. For example, if two noninbred individuals have two identicalbydescent alleles (S_{7}), either of two identityinstate patterns (
Recall that the set of parameters Δ= (Δ_{1}, Δ_{2},... Δ_{9}) correspond to the probabilities of each identitybydescent mode and completely quantify the degree of relatedness between individuals. Following Thompson (1975), the probability of observing a particular allelic pattern,
Parameter space: The maximumlikelihood estimate of the set of Δ is found by searching over the parameter space until a maximum is found. In general an algebraic solution is impossible; as a result numerical methods are used. The implementation used here is based on a translation of the simplex method (Presset al. 1992), a hillclimbing optimization technique, into a set of C++ classes (B. G. Milligan, unpublished data). Although it is possible for such methods to identify local rather than global maxima in the likelihood function, plots of the likelihood surface and evaluation of the algorithm from distinct initial conditions suggest that multiple optima are unlikely to exist, in agreement with the analytical results of Thompson (1975).
A number of possibilities exist for defining the parameter space over which the optimization will be carried out. The complete parameter space is, of course, eightdimensional, corresponding to the nine distinct parameters Δ_{j} = Pr (S_{j}) constrained by the fact that they sum to unity. The immediate purpose, however, is to consider the case of a large noninbred population. In this instance only the last three parameters are nonzero. For the purposes of this analysis, the maximumlikelihood estimate is obtained by optimization within the twodimensional parameter space defined by the parameters (Δ_{7}, Δ_{8}, Δ_{9}) constrained by their sum being unity. It is meaningless to admit solutions outside this region as they correspond to undefined values for the probability of identitybydescent (Thompson 1975).
One of the useful features of maximumlikelihood estimators is that they can be readily adapted to a variety of situations. In some cases it may be known that individuals are either fullsibs or unrelated (or any other pair of degrees of relatedness; Mousseauet al. 1998). In such a case, the likelihood could be maximized within the onedimensional parameter space representing a continuum of linear combinations of those two degrees of relatedness. Alternatively, it may be known that individuals are either fullsibs, halfsibs, or unrelated or fullsibs, parentoffspring, or unrelated (Thomaset al. 2002). In such cases, the optimization can be done within the parameter space defined by appropriate linear combinations of these three degrees of relatedness. Thus, the basiclikelihood estimator described here can easily be adapted to a variety of different population structures simply by appropriate choice of parameter space.
Of primary concern is the statistical behavior of the likelihood estimator, especially in contrast to existing alternatives (Queller and Goodnight 1989; Liet al. 1993; Ritland 1996a; Lynch and Ritland 1999; Wang 2002). Because of the reliance on numerical optimization, the statistical behavior must be determined by simulation. The basic simulation process involved generating replicate genetic data sets for a pair of individuals under conditions of known relatedness, a specified number of sampled loci, and a known distribution of allelic variation at each locus. Analytical results based on appendix a were used to verify the simulations.
Methodofmoments estimators: The performance of 6 estimators of relatedness, quantified as the coancestry coefficient θ, was investigated. The likelihood estimator calculated θ from Equation 1 as the maximumlikelihood estimate of the identitybydescent probabilities, Δ. Five additional nonlikelihood estimators were considered as being representative of the diversity of the ones available. They represent 5 different means of using the similarity in allelic states between individuals to construct estimates of relatedness, and the ones tested by Van de Casteele et al. (2001) performed favorably under some conditions in their evaluation of 10 different estimators. These have all been described and compared previously (Lynch and Ritland 1999; Van de Casteeleet al. 2001; Wang 2002); consequently, rather than repeating their formulations here the reader is referred to the specific equations presented by Ritland (1996a), Lynch and Ritland (1999), and Wang (2002).
The most commonly used estimator is one published by Queller and Goodnight (1989), of which a number of variants are possible. The variant chosen here was the symmetric one obtained by averaging (rˆ_{xy} + rˆ_{yx})/2 (Lynch and Ritland 1999, Equation 11) across loci. A second estimator of relatedness (Liet al. 1993) is based directly on the pattern of shared alleles between two individuals and was obtained by averaging rˆ_{xy} (Lynch and Ritland 1999, Equation 8) across loci. The third estimator, a methodofmoments one based on the correlation between relatedness and genotypic similarity (Ritland 1996a), was obtained as
Several of these estimators have undesirable behavior under certain conditions. For example, with two alleles the Queller and Goodnight (1989) estimator is undefined for heterozygous reference individuals, and for two equally frequent alleles the Lynch and Ritland (1999) estimator is also undefined for heterozygous reference individuals. Consequently, some loci must be discarded for multilocus estimates under these conditions, depending on the sampling of alleles at each locus. Both of these estimators are also based on arithmetic averages using each of the two individuals as a reference. At a particular locus one, both, or neither of those individuals can be heterozygous. Thus, the locusspecific value actually used in averaging across loci was the average (as described above) if both locusspecific values are defined and the defined value if only one is defined; otherwise, the locus was ignored. This approach attempts to maximize the amount of information obtained from these estimators within the constraints imposed by their mathematical definition.
Additionally, the methodofmoments estimators are generally not constrained to lie within the biologically relevant range of [0, 0.5], unlike the traditional maximumlikelihood estimator. This property enables them to remain statistically unbiased; however, individual estimates may not have meaning when interpreted as probabilities of identitybydescent. One means of handling this is to truncate the methodofmoments estimates to lie within the proper range, that is, to replace lower values with zero and larger values with 0.5. To investigate the effect of the parameter range itself, as opposed to the type of estimator, all methodofmoments estimators were examined in both the standard and the truncated form.
Simulations: A range of sampling conditions was considered to mimic the variety of different genetic markers that might be available for estimating relatedness. Although for some organisms huge arrays of polymorphic genetic loci are available, for the vast majority of natural populations this is not the case. Thus, the range in number of loci (530) mimics moderate genetic samples. A great variety of types of genetic markers are available for quantifying relatedness. These range from markers that segregate few distinguishable alleles, generally including allozymes, singlenucleotide polymorphisms, or restriction/PCR fragments, to markers that segregate many distinguishable alleles, commonly microsatellites. To reflect this range in marker types, loci segregating 2, 5, 10, and 20 alleles were considered. Three different allelefrequency distributions were used for the simulations: one in which all alleles occur at equal frequency, one in which a single allele occurs with a frequency of 0.8 and the remainder are equally frequent, and one in which allele frequencies at each locus were independently drawn from the same Dirichlet distribution (Stuart and Ord 1987, Exercise 5.33, p. 209) with all parameters set to unity. The last case approximates natural situations better by allowing the distributions to vary across loci; however, the first two are useful for isolating the effects of each factor. Finally, four actual degrees of relatedness between individuals were considered: parentoffspring, fullsibs, first cousins, and unrelated individuals. This range of conditions in relatedness, numbers of loci, numbers of alleles, and types of allelefrequency distributions was chosen to make the diversity of natural situations tractable, so that the influence of each main characteristic of a genetic sample on the statistical performance of the estimators can be examined. Future studies can focus more specifically on particular types of genetic markers or population structures.
To determine the statistical behavior of each estimator under each condition, sets of 1000 replicate samples of two individuals were obtained. Each of the six estimators was used to estimate θ for each of the replicate samples. The mean and standard error of the population of estimates were calculated from these samples. The bias of each estimator was quantified as the deviation of the mean from the known parametric value of θ used to generate the data. The root meansquare error was quantified as
RESULTS
As with any estimator of genetic relatedness, the quality of the estimate depends on the amount of available genetic information. Typically, both the number of loci for which genetic information is available and the number of alleles segregating at those loci have strong influences on the standard error of the estimate of relatedness. Additionally, different estimators of relatedness often respond differently to the amount of genetic information available.
Figure 2 illustrates the general level of variation yielded by each of the estimators. It is evident that the likelihood estimator described here has lower standard error than any of the others under all conditions; however, the other five are rather similar overall. One feature of the other estimators is their propensity to yield estimates of θ that lie outside the biologically meaningful range of [0, 0.5]. Under some conditions approximately half of the estimates are negative, for example. Because the focus of interest for these measures is on a specific pair of individuals, it is difficult to interpret the meaning of estimates that suggest, for example, that two individuals are less related than unrelated. However, because the likelihood estimator is constrained to always produce estimates within the biologically meaningful range, some bias is introduced near the boundary. For example, for unrelated individuals θ= 0, yet the likelihood estimator evidently commonly generates values that overestimate that. Thus, while exhibiting less variation, the likelihood estimator is more biased under some conditions. Clearly, the truncated estimators will also exhibit less variation and more bias than the untruncated ones for the same reason.
An additional feature that is evident from Figure 2 is that many of the estimators are skewed. This is particularly the case for the Queller and Goodnight (1989) and especially the Ritland (1996a) estimators, which interestingly are skewed in opposite directions. This skew may have a significant impact on the use of these estimators, because even though they are essentially unbiased in expectation, the modal estimate does not equal the expectation; the most probable outcome in any particular case will be an incorrect estimate.
Standard error: Figure 3 quantifies the standard error of each estimator of θ as a function of the amount of genetic information available. The variation for all estimators declines with the number of loci sampled. Generally, the standard error of the likelihood estimator is lower than that of any of the others. Interestingly, the estimators proposed by Ritland (1996a) and Lynch and Ritland (1999) approach the likelihood estimator under conditions of both low degree of relatedness and very restricted genetic information. Undoubtedly, this is because the weights used across loci for these estimators assume no relatedness. However, their performance is not consistent even for unrelated individuals; when more genetic information is available in the form of more uniform allelefrequency distributions, their performance is consistent with the other nonlikelihood estimators. The Ritland (1996a) estimator even performs distinctly worse than any other under some conditions considered. This anomalous behavior is because it involves a ratio of the allelic similarity (a number from the set {0, 0.25, 0.5, 1}) and the allele frequency. If by chance a rare allele is shared between individuals, this ratio can be substantial, greatly inflating the variance of the estimator and causing the extreme skewness observed in Figure 2. Thus, the utility of this estimator is highly dependent on the unknown quantity being estimated, relatedness, and the details of the genetic sample. To a lesser degree, the same is true of the Lynch and Ritland (1999) estimator. In contrast, the likelihood estimator maintains a low standard error across the full range of conditions. As such, it may be preferable overall, because in general no information is available to enable one to choose among the estimators on the basis of their performance under the special conditions of actual relatedness applying to a particular pair of individuals.
In these simulations there is no indication that the two different versions of the Wang (2002) estimators differ substantially. Both are essentially unbiased under all conditions and only for parentoffspring pairs were the differences between them >1%. Nor is there a strong indication that they represent substantial improvements over the other estimators. In no case is either better than the likelihood estimator. Perhaps these results are a consequence of the slightly different statistics used to evaluate the estimators. Whereas Lynch and Ritland (1999) and Wang (2002) use the singlelocus variance, quantified as a mean for 10locus samples, to evaluate the estimators, Figure 3 evaluates them on the basis of the actual variance for many different sample sizes. This approach is adopted for two important reasons. First, it more closely approximates the information needed directly for experimental design purposes, when the fundamental decision is often based on the reduction of variance as a function of genetic sample size. Second, the relative performance of the different estimators is a function of the number of loci sampled; basing an evaluation on a single sample size may be misleading for extrapolations to another.
The actual degree of relatedness between individuals has relatively little effect on the standard error of the likelihood estimator of θ. For example, the standard errors of 30locus likelihood estimates for fullsibs and first cousins differ by <4%, despite these representing quite different degrees of relatedness. The standard error is lower for unrelated individuals because of the constraint that estimates must be within the biologically realistic range. This independence of actual relatedness for standard error of θ is broadly consistent across all the estimators when the allelefrequency distribution is favorable; it is somewhat less so when the allelefrequency distribution is dominated by a single allele segregating at high frequency or when variation among loci in allelefrequency distributions exists.
Bias: The second main statistical feature of each estimator, bias, is illustrated in Figure 4. Two features are immediately apparent. All of the nonlikelihood estimators are essentially unbiased under all conditions; in contrast, the likelihood estimator is biased under some conditions. As mentioned before, this is a consequence of the fact that the likelihood estimator is constrained within the biologically meaningful range of [0, 0.5].
Unlike for standard error, the actual degree of relatedness does influence the bias of the likelihood estimator. For example, the likelihood estimator for parentoffspring and fullsib relationships yields estimates that are quite close to the true value of θ; in fact, its bias is either zero or close enough as to be biologically insignificant. However, when actual relatedness is close to the boundary, which is the case for both first cousins and unrelated individuals, the bias is much larger. This is true for any of the allelefrequency distributions. When the allelefrequency distribution is favorable, increasing numbers of loci can substantially reduce the bias of the likelihood estimator, potentially to the point of being biologically insignificant even for unrelated individuals. In contrast, when the allelefrequency distribution is dominated by a single allele, increasing numbers of loci have little effect on bias for reasonable numbers of loci.
Root meansquare error: The third main statistical feature of each estimator, root meansquare error, is illustrated in Figure 5. This measure is a reflection of the mean deviation of the distribution of estimates from the parametric value of θ used in the simulation. As such it integrates both the standard error and the bias of the estimators. Largely these curves follow the corresponding ones for standard error, an indication that under most conditions all estimators are essentially unbiased. Only in situations both lacking in useful genetic information (e.g., a single predominant allele at each locus) and with true relationships near the boundary of the parameter space (e.g., especially unrelated individuals) does the likelihood estimator perform notably worse than the others with regard to root meansquare error. In all other cases considered, the likelihood estimator performs better than any alternative with regard to this integrated measure of performance.
Truncated estimators: The performance of the truncated methodofmoments estimators relative to the likelihood estimator is largely anticipated from Figure 2. In the cases of fullsib and parentoffspring relatedness relatively few estimates lie beyond the meaningful range of [0, 0.5], so truncation has little effect; in contrast, for first cousins and unrelated individuals, substantial numbers of estimates do so and truncation has a large effect. For example, although the standard error is essentially unchanged between truncated and nontruncated estimators for the former two, it is somewhat reduced for first cousins and substantially reduced for unrelated individuals. For first cousins the truncated methodofmoments estimators exhibit no standard error lower than that of the maximumlikelihood estimator, and in most cases still exceed it. The reduction in standard error is great enough for unrelated individuals that, for the Dirichlet distribution of allele frequencies and when one allele predominates the Ritland (1996a) and Lynch and Ritland (1999) estimators, both exhibit a lower standard error than that of the maximumlikelihood estimator. The others, and all truncated estimators when alleles are equally frequent, exhibit standard errors similar to the maximumlikelihood estimator. Interestingly, the relative ranking of the methodofmoments estimators is quite similar whether they are truncated or not.
The bias of the truncated methodofmoments estimators increases substantially for the two lessrelated cases. Although the relative increase in bias is quite large (>10 to 100fold in some cases), the bias is somewhat less than that exhibited by the maximumlikelihood estimator. Under these conditions the truncated Ritland (1996a) and Lynch and Ritland (1999) estimators exhibited the lowest bias. However, the truncated Ritland (1996a) estimator is notably more biased under the two closer degrees of true relatedness, even exceeding the bias of the maximumlikelihood estimator under some conditions. This is another manifestation of the sensitivity of this estimator to the conditions under study, at least some of which will be unknown in any realistic situation.
Segregating alleles: The previous discussion illustrated the performance of alternative estimators under several different sampling conditions. These results are quite typical. However, the difference between the likelihood estimator and the others declines as the number of alleles segregating at each locus increases. That is, the nonlikelihood estimators improve in performance and approach the likelihood estimator as the amount of genetic information increases.
The variation of the likelihood estimator itself also changes in response to increasing numbers of alleles segregating at each locus (Figure 6). In all cases, more segregating alleles reduce the standard error of the likelihood estimate of θ. This is especially true when alleles are equally frequent. When one allele predominates, even a few additional alleles substantially reduce the standard error, which is not further reduced by many additional alleles. Thus, for a wide range of conditions a large reduction in standard error relative to biallelic samples is possible by sampling loci with even a few alleles; additional reductions are possible only if many alleles segregate at intermediate frequencies.
In contrast, the bias of the likelihood estimator is relatively unaffected by the number of alleles segregating at each locus (Figure 7). In fact, when the allelefrequency distribution is dominated by a single allele, the number of additional alleles segregating (and, as noted above, the number of loci) has essentially no influence on the bias over the range of loci considered (although the estimator is asymptotically unbiased). These sampling conditions are basically uninformative, so large samples are unhelpful. However, if the allelefrequency distribution is more favorable, the number of alleles segregating at each locus does have an influence on the bias of the likelihood estimator. As with the standard error, substantial reductions in bias are possible under all conditions of actual relatedness when even a few alleles are segregating. Even for conditions exhibiting the largest bias (e.g., unrelated individuals), the degree of bias can be reduced to biologically insignificant levels for realistic samples.
DISCUSSION
Despite its basic importance for understanding the biology of natural populations, estimation of relatedness between individuals remains a difficult challenge. A diversity of estimators (Thompson 1975; Queller and Goodnight 1989; Liet al. 1993; Ritland 1996a; Lynch and Ritland 1999; Wang 2002) have been developed to use the information contained within samples of molecular markers to estimate relatedness. With one exception (Thompson 1975), none of the methods proposed to date are traditional maximumlikelihood estimators. Here we investigate its statistical properties in comparison with five of the commonly used alternatives.
The prominent feature of the likelihood estimator is that it exhibits a lower standard error than that of any of the others for a wide range of reasonable sampling conditions. This conclusion appears to be independent of the number of loci sampled, the number of alleles segregating at each locus, or the frequency distribution of the segregating alleles. When many loci are sampled, the other estimators approach the performance of the likelihood estimator. This is especially true for the Ritland (1996a) and Lynch and Ritland (1999) estimators when individuals are only distantly related and the allelefrequency distribution is highly skewed. However, especially the Ritland (1996a) estimator is much less consistent across different sampling conditions and exhibits a strong dependency on the actual degree of relatedness, the unknowable quantity that is being estimated. Furthermore, the standard error of the likelihood estimator is largely independent of the actual degree of relatedness between individuals. Thus, from this standpoint the likelihood estimator exhibits preferable statistical behavior compared with any of the alternatives, because it is the only one that consistently maintains a low standard error across all conditions.
These conclusions differ dramatically from those obtained earlier for different maximumlikelihood estimators of relatedness (Ritland 1996a; Lynch and Ritland 1999). Their estimators performed so poorly as to be immediately discarded as useless in practice. This difference arises either from admitting solutions that cannot be directly interpreted biologically as probabilities of identitybydescent or from the nature of the likelihood function (see appendix b). In contrast, the likelihood estimator investigated here is consistent with the traditional literature on likelihood estimation of relatedness (Thompson 1975) and admits only solutions that are fully interpretable biologically. As a result, it performs much better than previously suggested for maximumlikelihood estimators.
The other feature of the likelihood estimator is that it, unlike the others, is biased under some conditions. The degree of bias is dependent on the actual degree of relatedness between individuals and the nature of the genetic information. If the actual degree of relatedness is near a boundary, such as for first cousins or unrelated individuals, the bias is more severe than if the actual degree of relatedness is within the interior of the parameter space, such as for fullsibs. However, the degree of bias can be greatly reduced by sampling loci that segregate for more alleles. Additional segregating alleles are not helpful if their frequency distribution is highly skewed. Samples of 2030 microsatellite loci, which often segregate for 20 or more alleles and exhibit high heterozygosity, could yield estimates of θ that are quite unbiased, even for unrelated individuals. However, even markers segregating for only a few alleles can also dramatically reduce the bias compared with biallelic markers. Thus, even though the likelihood estimator is more biased than any of the others, suitable genetic sampling can greatly reduce this problem. The amount of genetic information required to reduce the bias of the likelihood estimator to insignificant levels is within the range of feasible sampling efforts.
These two quantities, standard error and bias, are integrated by the measure of root meansquare error. As with the standard error, the likelihood estimator maintains a low root meansquare error under almost all conditions. The exceptions are cases involving little useful genetic information and true relationships that are near the boundary of the parameter space. However, the fact that a low root meansquare error is maintained by the likelihood estimator, despite the inherent potential for enhanced bias, indicates that in fact the bias is of little biological consequence. The other estimators more than make up for generally lower bias through their greater standard error. Additionally, the skewness identified for several other estimators may lead to further problems in practice.
Often the primary interest lies in the estimate of relatedness itself. Such would be the case, for example, in ascertaining family membership or determining the spatial structure of relatedness for sessile organisms. In such cases it is critical that each estimate of relatedness yield a biologically meaningful value. Such is the case for the likelihood estimator, which is constrained to yield estimates of θ within the biologically meaningful range of [0, 0.5]. The methodofmoments estimators may also be truncated to lie within the same range. Under conditions of low relatedness, this reduces their standard error and increases their bias. In many cases, however, the maximumlikelihood estimator still exhibits lower standard error or root meansquare error than that of the truncated ones. Thus, the general performance characteristics of this estimator are not strictly the result of constraining the parameter space to include only biologically meaningful estimates.
The primary interest in estimating relatedness may alternatively lie in using the estimates in a subsequent analysis. Such would be the case, for example, in estimating heritability or additive genetic variance from relatedness estimates and phenotypic observations (Ritland 1996b, 2000; Ritland and Ritland 1996; Mousseauet al. 1998; Thomas et al. 2000, 2002). In these cases, too, it may be problematic to allow relatedness estimates that lie beyond the biologically meaningful range. However, the amount of variation in estimates of θ can also be especially important, because of the propagation of error in the estimate of θ to variation in derived estimates of additive genetic variance, V_{A}.
Depending on the sampling conditions, standard errors for the nonlikelihood estimators are between 2 and 250% larger than the standard error for the likelihood estimator. That discrepancy is substantially improved only by truncating the nonlikelihood estimators when individuals are unrelated; unfortunately, in practice it will be unknown whether the specific pair of interest is unrelated or not.
This difference in standard error between estimators is likely to be significant. For example, in the study by Ritland and Ritland (1996) involving a sample of eight loci in Mimulus, the actual variance of relatedness was estimated to be only 0.04; almost all the observed variance in estimates of r was due to sampling error. Indeed, the general lack of application of markerbased estimates of relatedness to subsequent estimation of quantitative genetics parameters such as heritability or additive genetic variance may be attributable to the relatively poor performance of the available estimators. An estimator of relatedness, like the likelihood one discussed here, that exhibits both lower standard error and lower root meansquare error across many sampling conditions would improve the discriminatory power of such analyses.
Overall, the likelihood estimator discussed in this article offers several advantages compared with existing, commonly used estimators. First, except for some truncated estimators applied to unrelated individuals, it uniformly exhibits lower variation, even under conditions of relatively abundant genetic information. For example, even when 30 loci, each segregating for 20 equally frequent alleles, are sampled, the standard error is 10% (and >250% under some conditions) greater for some estimators than for the likelihood estimator. Second, all likelihood estimates are constrained to lie within the biologically meaningful range. Thus, the biological interpretation of individual estimates is quite straightforward. Although apparently not general practice, this constraint can be obtained by truncating the nonlikelihood estimates. While each estimate is now interpretable biologically, the statistical behavior of the truncated estimators is not generally improved over the maximumlikelihood estimator and in some cases is worse. Finally, the likelihood estimator naturally accommodates different genetic sampling schemes. For example, the relative weighting of data from microsatellite loci segregrating for many alleles in contrast to data from singlenucleotide polymorphisms segregating for only two alleles is accomplished directly by the likelihood function. Consequently, all available data can be used to ascertain the degree of relatedness between two individuals.
The main drawback associated with the likelihood estimator is that it can be biased under some conditions, especially if the true degree of relatedness is near the boundary of the parameter space and little genetic information is available. The same is also true of the truncated estimators, some of which are more biased than the likelihood function even when the true relatedness is not near the boundary. Even though biased, however, the maximumlikelihood estimator exhibits a lower root meansquare error than do alternative estimators under many conditions. Thus, the bias is quite minor from a biological perspective. Furthermore, the extent of the bias can be greatly reduced by suitable genetic sampling. If markers with many alleles segregating at intermediate frequencies are used, the bias can be reduced considerably. Often markers with only a few alleles segregating at intermediate frequencies approach the performance of highly polymorphic ones. Given the relative ease with which genetic information can be obtained, bias is not likely to be a major drawback in practice. However, a preliminary study aimed at defining the allelefrequency distributions prior to selecting loci for more intense sampling can dramatically reduce the genetic information required to obtain estimates of relatedness.
Although the maximumlikelihood estimator of relatedness performs extremely well overall, there are conditions under which others perform better according to specific metrics. None perform uniformly better under all conditions according to all metrics. Whether the most relevant performance metric is the standard error, the bias, the root meansquare error, or some other one likely depends in detail on the ultimate use of the relatedness estimates. For specific applications under specific genetic conditions it may be possible to identify one estimator that is optimal. However, for its wide range of applicability and excellent performance across almost all conditions it is difficult to improve on the maximumlikelihood estimator.
APPENDIX
A Although the likelihood function is in general complex and does not admit full analysis, some insight can be obtained for special cases and those results can be used to test the simulations. The first special case to consider corresponds to a single locus segregating for n equally frequent alleles. Table A1 gives a summary of this situation. Each of the allelic identityinstate patterns (see Table 1)
For the case of a single locus segregating for n equally frequent alleles, the maximum of the likelihood function can be solved analytically. In some cases (e.g.,
Even from Table A1 it is possible to understand the behavior of the maximumlikelihood estimator under more general conditions. For example, if an infinite number of alleles segregate in a noninbred population the identityinstate pattern
The effect of fewer alleles segregating is also evident from Table A1. In this case identityinstate patterns
The rapid improvement in performance of the maximumlikelihood estimator with increasing number of segregating alleles is also understandable from Table A1. With only two alleles, the set of observable identityinstate patterns is quite constrained and the maximumlikelihood estimates are less concordant with the actual mode of identitybydescent. Even one additional allele greatly improves the concordance.
A second special case that is amenable to analysis corresponds to a parentoffspring pair assayed at L loci, each of which segregates for n equally frequent alleles. In this case the true relationship is S_{8} = 1 and the likelihood function is given by
APPENDIX B
The methodofmoment estimator used by Lynch and Ritland (1999) is apparently from the same probability model as described above, based on both its derivation and its performance in simulations. However, the proposedlikelihood function (their Equation 12) differs from that presented by Thompson (1975) and generalized here in Table 1. The following illustrates this incompatibility beginning with the compact notation used by Ritland (2000), which expresses as a single expression the likelihood function for samples from a noninbred population (Thompson 1975).
First, let
Acknowledgments
Thanks are due to Kelly Gallagher, Kermit Ritland, Elizabeth Thompson, and two anonymous reviewers for their helpful comments. Thanks are also due to Jinliang Wang for his comments and for clarifying the details of the weighting scheme used in his estimator. Finally, Sara Knott and Peter Visscher provided many stimulating discussions and insightful comments that greatly improved the manuscript. This work was supported by National Science Foundation grant DEB9904026.
Footnotes

Communicating editor: J. B. Walsh
 Received July 17, 2002.
 Accepted November 26, 2002.
 Copyright © 2003 by the Genetics Society of America