Rank-invariant estimation of inbreeding coefficients

The two alleles an individual carries at a locus are identical by descent (ibd) if they have descended from a single ancestral allele in a reference population, and the probability of such identity is the inbreeding coefficient of the individual. Inbreeding coefficients can be predicted from pedigrees with founders constituting the reference population, but estimation from genetic data is not possible without data from the reference population. Most inbreeding estimators that make explicit use of sample allele frequencies as estimates of allele probabilities in the reference population are confounded by average kinships with other individuals. This means that the ranking of those estimates depends on the scope of the study sample and we show the variation in rankings for common estimators applied to different subdivisions of 1000 Genomes data. Allele-sharing estimators of within-population inbreeding relative to average kinship in a study sample, however, do have invariant rankings across all studies including those individuals. They are unbiased with a large number of SNPs. We discuss how allele sharing estimates are the relevant quantities for a range of empirical applications.


INTRODUCTION
Allelic dependence at a locus is usually quantified by inbreeding coefficients for individuals or populations, with these measures referring either to correlations of allelic state indicators (Wright, 1922) or to probabilities of identity by descent, ibd, (Malécot, 1948). Here we use ibd and we have advocated allele-sharing estimators ( (Weir & Goudet, 2017), WG17 henceforth; (Goudet et al., 2018)) that are unbiased for individual and population inbreeding coefficients relative to average kinships among specified pairs of individuals. Estimators such as those in PLINK ( (Purcell et al., 2007) and GCTA (Yang et al., 2011), that use sample allele frequencies, confound inbreeding estimates by the averages of individual kinships. Our work recognizes the need to estimate inbreeding coefficients from many millions of SNP genotypes where likelihood methods may not be feasible and we employ moment-based methods.
There have been many published accounts of inbreeding estimation, including the recent evaluation of several methods by Alemu et al. (2021). Among those that refer to allele sharing, Li & Horvitz (1953) discussed an inbreeding estimator based on observed homozygosity, i.e., within-individual sharing of maternal and paternal alleles. They compared observed sharing to the value expected without inbreeding. They also constructed an estimator from the proportions of each allele type that were homozygous in a sample and gave an expression that was investigated further by Ritland (1996). Ritland used allele sharing within and between individuals and his inbreeding estimates assumed "independence or nearindependence" of individuals. If individuals are not independent, the rankings of his inbreeding coefficient estimates change with the sample. In WG17 we estimated inbreeding coefficients by comparing within-individual allele-sharing to average sharing between pairs of individuals in a sample. By not making explicit use of sample allele frequencies, we preserved the ranking of estimates across different samples and this is our central theme here. Ritland's individual-level inbreeding coefficients were also derived by Yang et al. (2011) as the correlation between uniting gametes and were expressed in terms of allele dosages for an individual and sample allele frequencies. This estimator was written asF UNI in Yengo et al. (2017), and is less biased than the estimator in Yang et al. (2011) obtained from the diagonal elements of a genomic relationship matrix (GRM) of VanRaden (2008). We compare these two estimates below with allele-sharing and other methods: pedigree-based path-counting (Wright, 1922), maximum-likelihood estimation, MLE, (e.g., (Hall et al., 2012)) and runs of homozygosity (ROH) (e.g., (Ceballos et al., 2018)).

METHODS Statistical sampling
We can describe the dependence between pairs of uniting alleles in a single population without invoking an evolutionary model for the history of the population. In this "statistical sampling" framework (Weir, 1996) we do not consider the variation associated with evolutionary processes but we do consider the variation among samples from the same population. Although extensive sets of genetic data allow individuallevel inbreeding coefficients to be estimated with high precision, we start with population-level estimation.
Allelic dependencies can be quantified with the within-population inbreeding coefficient, written here as f W to emphasize it is a withinpopulation quantity, defined by where H l is the population proportion of heterozygotes for the reference allele at SNP l and p l is the population proportion of that allele. The same value of f W is assumed to apply for all SNPs. An immediate consequence of this definition is that the population proportions of homozygotes for the reference and alternative alleles are p 2 l þ p l ð1 À p l Þf W and ð1 À p l Þ 2 þ p l ð1 À p l Þf W respectively. This formulation allows f W to be negative, with the maximum of −p l /(1 − p l ) and −(1 − p l )/p l as lower bound. It is bounded above by 1. Hardy-Weinberg equilibrium, HWE, corresponds to f W = 0 and textbooks (e.g., (Hedrick, 2000)) point out that negative values of f W indicate more heterozygotes than expected under HWE.
Observed heterozygote proportionsH l have H l as within-population expectation E W over samples from the study population, E W ðH l Þ ¼ H l , and this would provide a simple estimator of f W if the population allele proportions were known. In practice, however, these proportions are unknown. Steele et al. (2014) suggested use of data external to the study sample to provide reference allele proportions in forensic applications where a reference database is used for making inferences about the population relevant for a particular crime. The more usual approach is to use study sample proportionsp l in place of the true proportions p l , as in equation 1 of Li & Horvitz (1953): The moment estimator in Eq. (2) is also an MLE of f W when only one locus is considered, but it is biased (Robertson & Hill, 1984) since not only is it a ratio of statistics but also the expected value E W ½2p l ð1 Àp l Þ over repeated samples of n from the population is 2p (e.g., (Weir, 1996), p39). This approach can be used to estimate the within-population inbreeding coefficient f j for each individual j in a sample from one population. These are the "simple" estimators of Hall et al. (2012) and thê f HOMj of Yengo et al. (2017): The sample heterozygosity indicatorH jl is one if individual j is heterozygous at SNP l and is zero otherwise. Averaging Eq.
(3) over individuals gives the estimator based on SNP l in Eq. (2).
A single SNP provides estimates that are either 1 or a negative value depending onp l , so many SNPs are used in practice. In both Hall et al. (2012) and Yengo et al. (2017) data were combined over loci as weighted or "ratio of averages" estimators: Gazal et al. (2014) referred to this estimator as f PLINK as it is an option in PLINK. We show below the good performance of this weighted estimator for large sample sizes and large numbers of loci. We will consider throughout that a large number L of SNPs are used so that ratios of sums of statistics over loci, such as in Eq. (4), have expected values equal to the ratio of expected values of their numerators and denominators. Ochoa & Storey (2021) showed statistics of the formÃ L =B L , whereÃ L ¼ P L l¼1 a l =L andB L ¼ P L l¼1 b l =L, have expected values that converge almost surely to the ratio A/B when E W ðÃ L Þ ¼ Ac L and E W ðB L Þ ¼ Bc L . This result rests on the expectations E W ða l Þ ¼ Ac l and E W ðb l Þ ¼ Bc l with c L ¼ P L l¼1 c l =L. It requires |a l |, |b l | to both be no greater than some finite quantity C, c L to converge to a finite value c as L increases, and for Bc not to be zero. For the ratio in Eq. (4), a l ¼H jl , b l ¼ 2p l ð1 Àp l Þ so A = (1 − f j ), B = 1 for large sample sizes n, and c L = ∑ l 2p l (1 − p l )/L ≤ 1/2. The conditions are satisfied providing at least one SNP is polymorphic. For an "average of ratios" estimator of the form P L l¼1 ða l =b l Þ=L, the denominators b l can be very small and convergence of its expected value is not assured.
As an alternative to using sample allele frequencies, Hall et al. (2012) used maximum likelihood to estimate population allele proportions for multiple loci whereas Ayres & Balding (1998) used Markov chain Monte Carlo methods in a Bayesian approach that integrated out the allele proportion parameters. Neither of those papers considered data of the size we now face in sequence-based studies of many organisms, and we doubt the computational effort to estimate, or integrate over, hundreds of millions of allele proportions in Eqs. (2) or (4) adds much value to inferences about f. The allele-sharing estimators we describe below regard allele probabilities as unknown nuisance parameters and we show how to avoid estimating them or assigning them values. Hall et al. (2012) used an EM algorithm to find MLEs for f j when population allele proportions were regarded as being known and equal to sample proportions. Alternatively, a grid search can be conducted over the range of validity for the single parameter f j that maximizes the log-likelihood Estimation of the within-population inbreeding coefficients f W (F IS of (Wright, 1922)) and f j does not require any information beyond genotype proportions in samples from a study population, nor does it make any assumptions about that population or the evolutionary forces that shaped the population. The coefficients are simply measures of dependence of pairs of alleles within individuals.

Genetic sampling
Inbreeding parameters of most interest in genetic studies are those that recognize the contribution of previous generations to inbreeding in the present study population. This requires accounting for "genetic sampling" (Weir, 1996) between generations, thereby leading to an ibd interpretation of inbreeding: ibd alleles descend from a single allele in a reference population. It also allows the prediction of inbreeding coefficients by path counting when pedigrees are known (Wright, 1922). If individual J is ancestral to both individuals j 0 and j″, and if there are n individuals in the pedigree path joining j 0 to j″ through J, then F j = ∑(0.5) n (1 + F J ) where F J is the inbreeding coefficient of ancestor J and F j is the inbreeding coefficient of offspring j of parents j 0 and j″. The sum is over all ancestors J and all paths joining j 0 to j″ through J. The expression is also the coancestry θ j 0 j 00 of j 0 and j″: the probability an allele drawn randomly from j 0 is ibd to an allele drawn randomly from j″.
The allele proportion p l in a study population has expectation π l over evolutionary replicates of the population from an ancestral reference population to the present time. Sample allele proportionsp l provide information about the population proportions p l , and their statistical sampling properties follow from the binomial distribution. We do not invoke a specific genetic sampling distribution for the p l about their expectations π l although we do assume the second moments of that distribution depend on probabilities of ibd for pairs of alleles. One consequence of the assumed moments is that the probability of individual j in the study sample being heterozygous, i.e., the total expected value E T of the heterozygosity indicator over replicates of the history of that individual, is The quantity F j is the individual-specific version of F IT of Wright (1922) and we can regard it as the probability the two alleles at any locus for individual j are ibd. There is an implicit assumption in Eq. (5) that the reference population needed to define ibd is infinite and in HWE: there is probability F j that j has homologous alleles with a single ancestral allele in that population and probability (1 − F j ) of j having homologous alleles with distinct ancestral alleles there. In the first place, the single ancestral allele has probability π of being the reference allele for that locus and the implicit assumption is that two ancestral alleles are both the reference type with probability π 2 . This does not mean there is an actual ancestral population with those properties, any more than use of E T means there are actual replicates of the history of any population or individual, and we note that Eq. (5) does not allow higher heterozygosity than predicted by HWE. Nonetheless, the concept of ibd allows theoretical constructions of great utility and we now present a framework for approaching empirical situations. Inbreeding, or ibd, implies a common ancestral origin for uniting alleles and statements about sample allele proportionsp l require consideration of possible ibd for other pairs of alleles in the sample. The total expectation of 2p l ð1 Àp l Þ over samples from the population and over evolutionary replicates of the study population is ( (Weir, 1996), p176) where F W is the parametric inbreeding coefficient averaged over sample members, F W ¼ P n j¼1 F j =n, and θ S is the average parametric coancestry in the sample, θ S ¼ P n j¼1 P j 0 ≠j θ jj 0 =½nðn À 1Þ. Equivalent expressions were given by McPeek et al. (2004) and DeGiorgio and Rosenberg (2009). We note the relationship f W = (F W − θ S )/(1 − θ S ) given by Wright (1922) and we showed in WG17 the equivalent expression For a large number of SNPs, the expectation of a ratio estimator of the type considered here is the ratio of expectations (Ochoa & Storey, 2021). Therefore, the total expectations of thef Homj , taking into account both Q.S. Zhang et al. statistical and genetic sampling, are For all sample sizes,f HOMj has an expected value less than the true value f j , with the bias being of the order of 1/n. The ranking of E T ðf HOMj Þ values, however, is the same as the ranking of the f j and, therefore, of the F j . For large sample sizes, Eq. (7) reduces to E T ðf HOMj Þ ¼ f j . Averaging over individuals shows that E T ðf HOM Þ ¼ f W : the population-level estimator in Eq.
(2) has total expectation of f W , not F W . A different outcome is found for thef UNIj estimator of Yengo et al. (2017) (i.e.,f III of Yang et al. (2011);f GCTA3 of (Gazal et al., 2014)). This estimator, with the weighted (w) ratio of averages over loci we recommend, as opposed to the unweighted (u) average of ratios over loci used in their papers, iŝ In this equation X jl is the reference allele dosage, the number of copies of the reference allele, at SNP l for individual j. It is equivalent to the estimator given by (Ritland (1996), eq. 5) and attributed by him to Li & Horvitz (1953). Ochoa & Storey (2021) showed thatf w UNIj has expectation, for a large number of SNPs and a large sample size, of where Ψ j is the average coancestry of individual j with other members of the study sample: have an average of θ S over members of the sample, so the average of the ψ j 's is zero and expected value of the average of thef w UNIj is f W , as is the case forf ASj below.
Equation (9) shows that thef w UNIj have expected values with the same ranking as the F j values only if every individual j in the sample has the same average kinship ψ j with other sample members.
Finally, we mention another common estimator described by VanRaden (2008), termed f GCTA1 by Gazal et al. (2014) and available from the GCTA software (Yang et al., 2011) with option --ibc. We referred to this as the "standard" estimator in WG17. The weighted version for multiple loci iŝ and it has the large-sample expectation of (f j − 4ψ j ) as is implied by WG17 (Eq. 13) and as was given by Ochoa & Storey (2021). We summarize the various measures of inbreeding and coancestry in Table 1, and we include sample sizes in the expectations shown in Table 2.
Thef HOM ,f UNI ;f STD andf MLE estimators of individual or population inbreeding coefficients make explicit use of sample allele proportions. This means that all four have small-sample biases, and none of the four provide estimates of the ibd quantities F or F j . We showed thatf HOM is actually estimating the within-population inbreeding coefficients: the total inbreeding coefficients relative to the average coancestry of pairs of individuals in the sample, butf UNI andf STD are estimating expressions that also involve average kinships ψ.

Allele sharing
In a genetic sampling framework, and with the ibd viewpoint, we consider within-individual allele sharing proportions A jl for SNP l in individual j (we wrote M rather than A in WG17 and in (Goudet et al., 2018)). These equal one for homozygotes and zero for heterozygotes and sample values can be expressed in terms of allele dosages,Ã jl ¼ ðX jl À 1Þ 2 . We also consider between-individual sharing proportions A jj 0 l for SNP l and individuals j and j 0 . These are equal to one for both individuals being the same homozygote, ibd probability for homologous alleles F Gold : Actual ibd in simulations. θ jj 0 Coancestry for individuals j; j 0 : ibd probability θ PED : Path counting.
for random alleles from j and j 0 . θ Gold : Actual ibd in simulations.
The following hold for PED and Gold values. No explicit expression.
f MLEj Maximization of likelihood for f j .
No explicit expression.
For weighted averages over large numbers of loci.
zero for different homozygotes, and 0.5 otherwise. Observed values can be written asÃ jj 0 l ¼ ½1 þ ðX jl À 1ÞðX j 0 l À 1Þ=2, with an average over all pairs of distinct individuals in a sample ofÃ Sl . Astle & Balding (2009) introducedÃ jj 0 l as a measure of identity in state of alleles chosen randomly from individuals j and j 0 , and Ochoa & Storey (2021) used a simple transformation of this quantity. The allele sharing for an individual with itself is A jjl = (1 + A jl )/2. The same logic that led to Eq. (5) provides total expectations for allelesharing proportions for all j; j 0 : Note that θ jj = (1 + F j )/2. The nuisance parameter 2π l (1 − π l ) cancels out of the ratio E T ðÃ jj 0 l ÀÃ Sl Þ=E T ð1 ÀÃ Sl Þ and this motivates definitions of allelesharing estimators of the inbreeding coefficient for individual j and the kinship coefficient for individuals j; j 0 aŝ For a large number of SNPs, these are unbiased for f j and ψ jj 0 for all sample sizes. We showed in WG17 there is no need to filter on minor allele frequency to preserve the lack of bias. Note thatf ASj is a linear function of the form a S þ b SÃj withÃ j being the total homozygosity for j and constants a S , b S being the same for all individuals j. Changing the scope of the study, from population to world for example, preserves linearity (with different values of a S , b S ). The changed estimates are linear functions of the old estimates: old and new estimates are completely correlated and are rank invariant over all samples that include particular individuals, i.e., over all reference populations. Unlike the case forf UNI orf STD , rank invariance is guaranteed forf ASj for any two individuals even if only one more individual is added to the study. For large sample sizes, ð1 ÀÃ Sl Þ % 2p l ð1 Àp l Þ. Under that approximation, f ASj is the same asf Homj but the approximation is not necessary in computerbased analyses. Summing the large-sample estimates over individuals not equal to j gives an estimator for the average individual kinship ψ j : Adding 2ψ ASj tof w UNIj givesf ASj , as expected, as does adding 4ψ ASj tô f w STDj . Similarly,ψ AS jj 0 is obtained by addingψ ASj andψ AS j 0 toψ STD jj 0 , where (Yang et al., 2011) ψ STD jj 0 ¼ P l ðX jl À 2p l ÞðX j 0 l À 2p l Þ P l 4p l ð1 Àp l Þ These are the elements of the first method for constructing the GRM given by VanRaden (2008). When inbreeding and coancestry coefficients are defined as ibd probabilities they are non-negative, but the within-population values f and ψ will be negative for individuals, or pairs of individuals, having smaller ibd allele probabilities than do pairs of individuals in the sample, on average. Individualspecific values of f always have the same ranking as the individual-specific F values, and they are estimable. Negative estimates can be avoided by the transformation to ðf ASj Àf min ASj Þ=ð1 Àf min ASj Þ wheref min ASj is the smallest value over individuals of thef ASj 's. We don't see the need for this transformation, and we noted above the recognition of the utility of negative values. Ochoa & Storey (2021) wished to estimate F j rather than f j and, to overcome the lack of information about the ancestral population serving as a reference point for ibd, they assumed the least related pair of individuals in a sample have a coancestry of zero. We showed in WG17 that this brings estimates in line with path-counting predicted values when founders are assumed to be not inbred and unrelated, but we prefer to avoid the assumption. We stress that, absent external information or assumptions, F is not estimable. Instead, linear functions of F that describe ibd of target pairs of alleles relative to ibd in a specified set of alleles are estimable and have utility in empirical studies.

Runs of homozygosity
Each of the inbreeding estimators considered so far has been constructed for individual SNPs and then combined over SNPs. Observed values of allelic state are used to make inferences about the unobserved state of identity by descent. Estimators based on ROH, however, suppose that ibd for a region of the genome can be observed. Although F is the probability an individual has ibd alleles at any single SNP, in fact ibd occurs in blocks within which there has been no recombination in the paths of descent from common ancestor to the individual's parents. Whereas a single SNP can be homozygous without the two alleles being ibd, if many adjacent SNPs are homozygous the most likely explanation is that they are in a block of ibd (Gibson et al., 2006). There can be exceptions, from mutation for example, and several publications give strategies for identifying runs of homozygotes for which ibd may be assumed (e.g., Gazal et al. (2014); (Joshi et al., 2015)). These strategies include adjusting the size of the blocks, the numbers of heterozygotes or missing values allowed per block, the minor allele frequency, and so on. These software parameters affect the size of the estimates (Meyermans et al., 2020). Some methods (e.g., Gazal et al. (2014); (Narasimhan et al., 2016)) use hidden Markov models where ibd is the hidden status of an observed homozygote. Model-based approaches necessarily have assumptions, such as HWE in the sampled population.
We provide more details elsewhere, but we note here that ROH methods offer a useful alternative to SNP-by-SNP methods even though they cannot completely compensate for lack of information on the ibd reference population. We note also that shorter runs of ibd result from more distant relatedness of an individual's parents, and ROH procedures can be set to distinguish recent (familial) ibd from distant (evolutionary) ibd. SNP-by-SNP estimators do not make a distinction between these two time scales.

Simulation study
We used the quantiNemo software (Neuenschwander et al., 2019) to simulate a five-generation pedigree of hermaphroditic individuals mating randomly, excluding selfing, with each mating producing a number of offspring drawn from a Poisson distribution with mean two. The zero-th generation was made of 50 founders, the first generation had 47 individuals and the second, third, fourth and fifth generations had 58, 56, 57, and 65 individuals respectively. This pedigree was then fed to a custom R script to draw gametes from each parent at each reproductive event, allowing for recombination based on a 20 Morgan recombination map with a genetic marker every 0.1 cM, for a total of 20,000 markers.
Each of the 100 alleles per marker among the 50 founders was given a unique identifier so that alleles in subsequent generations with the same identifier had actual identity by descent relative to the founders. The average actual ibd proportions over loci, within individuals and between each pair of individuals, provided "gold standard" inbreeding and coancestry coefficients, as opposed to the pedigree-based values we calculated by path counting. The gold values for inbreeding coefficients F j and coancestry coefficients θ jj 0 then allow calculation of gold values for f j , ψ j and, therefore, f STDj and f UNIj .
Finally, the two unique identifiers for each marker of the 50 founders were mapped to the SNP genotypes of the 50 founders generated with the msprime program (Kelleher et al., 2016) as follows: we assume the founders originated from a population with effective size N e = 10 4 , mutation rate μ = 10 −9 , recombination rate between neighboring base pairs r = 10 −7 . We assumed 20 chromosomes each 10 Megabase (10 7 ) long. The necessary arguments are mspms 100 20 -t 400 -r 40000 10000000 -p 9. This generated a dataset of 100 gametes and over 40,000 SNPs, with the first 20,000 used for the mapping of unique identifiers to SNP alleles. This mapping was applied to the genotypes of the non-founder individuals of the pedigree to generate their SNP genotypes.
The pedigree was constructed to provide fairly high levels of predicted coancestry among pairs of the 283 non-founder individuals, ranging from 0 to 0.464, with a mean of θ S = 0.053, assuming the 50 founders were unrelated and not inbred. The pedigree inbreeding coefficients ranged from 0 to 0.367, with a mean of F W = 0.050. The within-population inbreeding coefficient for the set of 283 non-founder individuals is f = (F W − θ S )/(1 − θ S ) = −0.003. Note, however, that the 50 individuals regarded as founders for the subsequent 283 had their own joint histories from the msprime simulation. These 50 had an average within-individual allele sharing ofÃ W ¼ 0:80385 and an average between-individual allele sharing ofÃ S ¼ 0:80355. The difference of these two proportions, Q.S. Zhang et al. which would be zero for a reference set of non-inbred and unrelated individuals, provides a within-founder allele-sharing inbreeding coefficientf W of 0.0015.
The various estimators of inbreeding examined with these data are shown in Table 2, and the correlation coefficients for each pair of estimates over the whole set of 283 non-founder individuals are shown in Table 3. There are very high correlations between pedigree and gold-standard values and also very high correlations betweenf HOM andf AS values, both as expected. There are lower correlations off UNI andf STD with pedigree-based or gold-standard inbreeding coefficients since those estimates reflect both f and ψ.
We see in Table 3 thatF ROH values are the most highly correlated with F Gold : this high correlation was obtained by adjusting the block size (100 SNPs) and the block overlap amount (50 SNPs) to bring estimates close to the known F Gold values. In practice the F Gold values are not known and the other estimators are all evaluated without external information. The high correlation off AS and maximum likelihood values suggests thatf MLE is estimating f rather than F because it uses the sample allele frequencies in place of the unknown allele probabilities. The weighted and unweighted versions off UNI are highly correlated with each other and with their gold values, but not with f Gold . There are generally low correlations for weighted and unweighted f STD values. Figure 1 (left) illustrates the linear relationship between f Pedj and F Pedj : f Pedj ¼ ðF Pedj À θ PedS Þ=ð1 À θ PedS Þ where θ PedS ¼ 0:053 is the average coancestry of pairs of non-founders, calculated from the pedigree. The F Goldj and f Goldj values are also correlated with the corresponding pedigree values, as is shown for f Goldj in Fig. 1 (center). The variation we see in Fig. 1 (center) for f Goldj around F Pedj reflects the variation of actual inbreeding about expected values, even for whole genomes, pointed out by Hill & Weir (2011). Wang (2016) showed that the number of SNPs also has an effect. The lack of relationship between pedigree-based values of individual average coancestry ψ j and individual inbreeding f j , leading to variable rankings for some estimators based on sample allele frequencies, is shown in Fig. 1 (right). Figure 2 (left) illustrates the similarity ofF ROH and F Gold and Fig. 2 (center) shows general agreement betweenF ROH andf AS , bearing in mind thatf AS estimates (F − θ S )/(1 − θ S ). Figure 2 (right) shows general agreement of the allele-sharing estimatorsf ASj with the gold-standard within-population inbreeding coefficients f Goldj . Figure 3 showsf UNIj to be a better estimator of f Goldj than isf STDj , as noted by Yang et al. (2011), and better performance for the weighted than unweighted averages over SNPs but still not as good asf ASj .

genomes data
We used 77m SNPs from the 22 autosomes for the 26 populations of the 1000 Genomes whole genome data to estimate inbreeding coefficients for all 2504 individuals in the project. Our focus was on the algebraic invariance of estimate rankings as the reference set of individuals changed from the population from which each individual was sampled, to the continental group for that population, to the whole world. We calculated the estimateŝ f ASj andf u UNIj for each individual and each reference set, and ranked estimates within each population. The two sets of   Figure 4 shows that within-population inbreeding coefficientsf AS for all 1000 Genomes populations outside the AMR group are essentially the same, and generally close to zero, when they are estimated relative to average coancestry within each population or continental group but change when the complete set of 26 populations is used as a reference. These latter values compare the allele sharing for each individual to the same reference value, the average sharing over all pairs of individuals in the whole dataset. The world reference gives markedly lowerf AS values for the African populations (AFR), reflecting their higher levels of genetic diversity. The rankings forf AS within a population, by construction, do not change with reference set. Highf AS values reflect admixture, consanguineous matings and high evolutionary coancestry. In contrast, thef UNI values are higher for African individuals than for any other individuals when the allele frequencies are from all 26 populations: this reflects an African-specific pattern of negative average individual kinships ψ, shown in the bottom row of Fig. 5.
The critical role that average kinship plays in inbreeding estimation is illustrated in Fig. 5. With each reference set, the allele-sharing inbreeding estimatesf AS are clustered for European (EUR) individuals, a little more diverse for East Asian (EAS) individuals, much more diverse for South Asian (SAS) and African (AFR) individuals, and extremely diverse for American (AMR) individuals. These values are consistent with those reported for the numbers of variant sites per genome (The 1000 Genomes Project Consortium, 2015). The variation among African and American average kinshipsψ AS is substantial: as these quantities determine how the expected values off UNI andf STD differ from the f target parameters, it is clear that these estimates cannot be used to rank individuals by their inbreeding levels.
For the African population ASW, individual NA20294 hasf AS values of −0.009, 0.001,−0.130 using ASW, AFR or World as a reference set and each estimate is ranked as number 16 among the 61 ASW estimates. The same individual hasf u UNI values of −0.007 (rank 36), 0.001 (rank 16) and 0.028 (rank 60) using ASW, AFR or World allele frequencies. Estimatorf u UNI indicates NA20294 to be among the least inbred of the ASW individuals when AFR sample allele frequencies are used, but among the most inbred when world-wide sample allele frequencies are used, even though the individual's own genotype is the same for each    analysis. Other examples of rankings changing with reference population forf UNI are shown in Fig. S3; for the admixed ACB and ACB populations, for example, the individuals appearing the most inbred with continental reference appear the least inbred with world reference and vice versa. This can have implications for studies of inbreeding depression, where trait values are regressed on estimated inbreeding coefficients.
A comparison of runs-of-homozygosity estimatesF ROHj with SNP-by-SNP estimates is shown in Fig. 6. The ROH estimates were produced with the --homozyg --homozyg-snp2 --homozyg-kb100 options in PLINK (Meyermans et al., 2020). The values of F ROHj depend on the PLINK settings for minor allele frequency pruning and linkage disequilibrium pruning, as well as on SNP density, so their expected values may differ from the true F j values. The left panel showsf ASj values and these have a correlation of 0.998 withF ROHj . The right panel showsf u UNIj estimates and these have a correlation of −0.337 withF ROHj estimates. Gazal et al. (2015) reported inbreeding estimatesF Fsuitej from ROH, although their method requires sample allele frequencies and so may have estimates of F confounded by average individual-specific average kinships. They also assumed Hardy-Weinberg equilibrium. However, there is good agreement off ASj values withF Fsuitej values (Fig. S4). The agreement between F Fsuitej andf u UNIj is seen there to be not as good.

DISCUSSION
Discussions on the estimation of individual inbreeding coefficients generally refer to F, the probability an individual has pairs of homologous alleles that are identical by descent. Among the estimators we have considered here,F ROH addresses F by assuming that long runs of homozygous SNPs represent ibd regions. The ROH estimates, however, are conditional on the settings used to calculate the estimates, and actual ibd in short runs of homozygotes may be ignored, so the expected values of these estimators is not known. The Bayesian approach of Vogl et al. (2002) also addresses F but at the computational cost of estimating allele proportions in a reference population assumed to have zero inbreeding or relatedness. All the other estimators considered here are, instead, addressing the within-population inbreeding coefficient f that compares F values to ibd probabilities for pairs of individuals. There is no need to specify the reference population implicit in the definition of identity by descent. There is also no need to assume the particular individuals in a sample have an inbreeding coefficient of zero. For large numbers of SNPs, allele-sharing estimatorsf AS are unbiased for f for all sample sizes and have values for a set of individuals that have invariant ranks over studies that include that set. We show that most estimators using sample allele frequencies are estimating some combination of f and of individual-specific average kinships ψ with individuals in the study. Estimators with expectations depending on ψ do not have invariant rankings, as we showed with data from the 1000 Genomes project as the study scope varied from the population to the continent to the world.
Our ibd-based model rests on expectations of allele-sharing proportions satisfying expressions such as Eq. (5). There is no requirement for nonoverlapping generations, or homogeneous populations, for example. This generality is a consequence of not needing allele frequencies, whether these refer to a population or to an individual.
The role of ibd probabilities in theoretical population and quantitative genetic contexts is well known, but we suggest it is rank-invariant estimators for the within-population parameters f j that are of relevance for empirical studies and we offer the examples in the following sections.

Genotype probabilities
There is often a need to estimate genotype probabilities from observed allele proportions using formulations with allele probabilities and ibd probabilities F (e.g., (National Research Council, 1996) for forensic science). Following Eq. (6) we see that it is 2p l ð1 Àp l Þð1 À f j Þ rather than 2p l ð1 Àp l Þð1 À F j Þ that is unbiased for 2π l (1 − π l )(1 − F j ) if F j and f j are known.

Inbreeding depression
Inbreeding is known to affect, linearly, the expected value of quantitative traits, and studies of inbreeding depression often proceed by regressing trait means on inbreeding levels. In Yengo et al. (2017), we usedF ROH ,f HOM andf UNI as inbreeding estimates and Kardos et al. (2018) pointed out that we did not discuss the distinction between F and f. We responded  (Yengo et al., 2018) with reasons for not wishing to useF ROH and we could have pointed out the linear relationship between f j and F j and the high correlation we showed above betweenf ASj andF ROHj means that regressing on eitherF ROH orf AS should lead to similar results. In less-homogeneous populations than represented by the UK Biobank data (Allen et al., 2012) we used in Yengo et al. (2017), it would appear to be better to usef ASj thanf UNIj to avoid any effects of individual-specific average kinships on inbreeding estimates. The correlation of trait andf ASj values is invariant over reference populations. Alemu et al. (2021) pointed out thatf HOM (andf AS ), gives equal weights to all SNPs, whereasf u UNI gives greater weight to SNPs with rare alleles. Alemu et al. did not consider the role of individual average kinships in the bias off UNI .

Genetic relatedness matrix
Inbreeding is also known to affect, linearly, the additive component of genetic variance. For additive traits, the genetic variance for individual j is ð1 þ F j Þσ 2 A where σ 2 A is the additive variance for populations in Hardy-Weinberg equilibrium. Consequently, the expected value of the sample varianceṼ T of trait values over a sample of n individuals is (Speed et al., 2012) Here the trait is additive and the errors, with variance σ 2 e , are independent of genetic effects. The GRM G has trace trðGÞ and sum of off-diagonal elements Σ G . If the GRM elements are (1 + F j ) on the diagonal and 2θ jj 0 off the diagonal then the trace is n(1 + F W ) and the sum of off-diagonal elements is n(n − 1)θ S so the genetic component of V T is ð1 þ F W À 2θ S Þσ 2 A . If the GRM is replaced by a matrix with allele-sharing inbreeding and kinship estimates, this becomes ð1 þ f W Þσ 2 A , reflecting that it is the within-population estimated GRM that is used in practice. We show elsewhere that the same expected variance holds with GRMs constructed withf STD orf UNI .
In summary, we have shown that inbreeding measures of utility in empirical studies are "within-population" with the choice of population being at the discretion of the investigator. With allele-sharing inbreeding estimators, the population specifies the set of individuals whose pairwise coancestry is the reference against which inbreeding is measured. For estimators making explicit use of sample allele frequencies, it is the population that furnishes those frequencies, although then inbreeding estimates are confounded by individual-specific average kinships. We showed algebraically and empirically that allele-sharing estimators have invariant rankings across choice of population.

DATA AVAILABILITY
The simulated data are available in the online supplement. The 1000 Genomes data are available at ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/.