[$\ast$]Corresponding author: Department of Mathematical Sciences,
Universitetsparken 5, 2100 Copenhagen, Denmark. Email<EMAIL_ADDRESS>
# Evaluation of population structure inferred by principal component analysis
or the admixture model

Jan van Waaij  Department of Mathematical Science, University of Copenhagen,
2100 Copenhagen, Denmark  Current address: Department of Health Technology,
Danish Technical University, 2800 Kgs. Lyngby, Denmark  These authors
contributed equally to this work  Song Li  Department of Mathematical Science,
University of Copenhagen, 2100 Copenhagen, Denmark  These authors contributed
equally to this work  Genís Garcia-Erill  Department of Biology, University of
Copenhagen, 2100 Copenhagen, Denmark  Anders Albrechtsen  Department of
Biology, University of Copenhagen, 2100 Copenhagen, Denmark  Carsten Wiuf
Department of Mathematical Science, University of Copenhagen, 2100 Copenhagen,
Denmark

###### Abstract

Principal component analysis (PCA) is commonly used in genetics to infer and
visualize population structure and admixture between populations. PCA is often
interpreted in a way similar to inferred admixture proportions, where it is
assumed that individuals belong to one of several possible populations or are
admixed between these populations. We propose a new method to assess the
statistical fit of PCA (interpreted as a model spanned by the top principal
components) and to show that violations of the PCA assumptions affect the fit.
Our method uses the chosen top principal components to predict the genotypes.
By assessing the covariance (and the correlation) of the residuals (the
differences between observed and predicted genotypes), we are able to detect
violation of the model assumptions. Based on simulations and genome wide human
data we show that our assessment of fit can be used to guide the
interpretation of the data and to pinpoint individuals that are not well
represented by the chosen principal components. Our method works equally on
other similar models, such as the admixture model, where the mean of the data
is represented by linear matrix decomposition.

###### keywords:

PCA; residuals; population modelling; ancient DNA; statistical fit

††articletype: inv

xx xx, xxxx xx xx, xxxx

##  1 Introduction

Principal component analysis (PCA) and model-based clustering methods are
popular ways to disentangle the ancestral genetic history of individuals and
populations. One particular model, the admixture model (Pritchard et al.,
2000), has played a prominent role because of its simple structure and, in
some cases, easy interpretability. PCA is often seen as being model free but
as noted by Engelhardt and Stephens (2010), the two approaches are very
similar. The interpretation of the results of a PCA analysis is often based on
assumptions similar to those of the admixture model, such that admixed
individuals are linear combinations of the eigenvectors representing unadmixed
individuals. In this way, the admixed individuals lie in-between the unadmixed
individuals in a PCA plot. As shown for the admixture model, there are many
demographic histories that can lead to the same result (Lawson et al., 2018a)
and many demographic histories that violate the assumptions of the admixture
model (Garcia-Erill and Albrechtsen, 2020). As we will show, this is also the
case for PCA, since it has a similar underlying model (Engelhardt and
Stephens, 2010).

The admixture model states that the genetic material from each individual is
composed of contributions from $k$ distinct ancestral homogeneous populations.
However, this is often contested in real data analysis, where the ancestral
population structure might be much more complicated than that specified by the
admixture model. For example, the $k$ ancestral populations might be
heterogeneous themselves, the exact number of ancestral populations might be
difficult to assess due to many smaller contributing populations, or the
genetic composition of an individual might be the result of continuous
migration or recent backcrossing, which also violates the assumptions of the
admixture model. Furthermore, the admixture model assumes individuals are
unrelated, which naturally might not be the case. This paper is concerned with
assessing the fit of PCA building on the special relationship with the
admixture model (Engelhardt and Stephens, 2010). In particular, we are
interested in quantifying the model fit and assessing the validity of the
model at the level of the sample as well as at the level of the individual.
Using real and simulated data we show that the fit from a PCA analysis is
affected by violations of the admixture model.

We consider genotype data $G$ from $n$ individuals and $m$ SNPs, such that
$G_{si}\in\\{0,1,2\\}$ is the number of reference alleles for individual $i$
and SNP $s$. Typically, $G_{si}$ is assumed to be binomially distributed with
parameter $\Pi_{si}$, where $\Pi_{si}$ depends on the number of ancestral
populations, $k$, their admixture proportions and the ancestral population
allele frequencies. For clustering based analysis such as ADMIXTURE (Alexander
and Lange, 2011), $k$ is the number of clusters while in PCA, it is the $k-1$
top principal components. We give the specifics of the admixture model in the
next section and show its relationship to PCA in the Material and methods
section.

Several methods aim to estimate the best $k$ in some sense (Alexander and
Lange, 2011; Evanno et al., 2005; Pritchard et al., 2000; Raj et al., 2014;
Wang, 2019), but finding such $k$ does not imply the data fit the model
(Lawson et al., 2018b; Janes et al., 2017). In statistics, it is standard to
use residuals and distributional summaries of the residuals to assess model
fit (Box et al., 2005). The residual of an observation is defined as the
difference between the observed and the predicted value (estimated under some
model). Visual trends in the residuals (for example, differences between
populations) are indicative of model misfit, and large absolute values of the
residuals are indicative of outliers (for example due to experimental errors,
or kinship). If the model is correct, a histogram of the residuals is expected
to be mono-modal centered around zero (Box et al., 2005).

In our context, Garcia-Erill and Albrechtsen (2020) argue that trends in the
residual correlation matrix carries information about the underlying model and
might be used for visual model evaluation. A method is designed to assess
whether the correlation structure agrees with the proposed model, in
particular, whether it agrees with the proposed number of homogeneous
ancestral populations (Garcia-Erill and Albrechtsen, 2020). However, even in
the case the model is correctly specified, the residuals are in general
correlated (Box et al., 2005), and therefore, trends might be observed even if
the model is true, leading to incorrect model assessment. To adjust for this
correlation, a leave-one-out procedure, based on maximum likelihood estimation
of the admixture model parameters, is developed that removes the correlation
between residuals in the case the model is correct, but not if the model is
misspecified (Garcia-Erill and Albrechtsen, 2020). This approach could also be
applied to PCA, where expected genotypes could be calculated using
probabilistic PCA (Meisner et al., 2021). This leave-one-out procedure is,
however, computationally expensive.

To remedy the computational difficulties, we take a different approach to
investigate the correlation structure. We suggest two different ways of
calculating the correlation matrix of the residuals. The first is simply the
empirical correlation matrix of the residuals. The second might be considered
an estimated correlation matrix, based on a model. Both are simple to compute.
Under mild regularity assumptions, these two measures agree if the model is
correct and the number of SNPs is large. Hence, their difference is expected
to be close to zero, when the admixture model is not violated. If the
difference is considerably different from zero, then this is proof of model
misfit.

To explore the adequacy of the proposed method, we investigate different ways
to calculate the predicted values of the genotype (hence, the residuals),
using Principal Component Analysis (PCA) in different ways. However, we also
show that this approach can be used on estimated admixture proportions.
Specifically, we use 1) an uncommon but very useful PCA approach (here, named
PCA 1) based on unnormalized genotypes (Cabreros and Storey, 2019; Chen and
Storey, 2015), 2) PCA applied to mean centred data (PCA 2), see Patterson et
al. (2006), and 3) PCA applied to mean and variance normalised data (PCA 3)
(Patterson et al., 2006). All three approaches are computationally fast and do
not require separate estimation of ancestral allele frequencies and population
proportions, as in Garcia-Erill and Albrechtsen (2020). Hence, the computation
of the residuals are computationally inexpensive. Additionally, we show that
this approach can also be applied to output from, for example, the software
ADMIXTURE (Alexander et al., 2009) to estimate $\Pi_{si}$ for each $s$ and
$i$, and to calculate the residuals from these estimates. An overview of PCA
can be found in Jolliffe and Cadima (2016).

We demonstrate that our proposed method works well on simulated and real data,
when the predicted values (and the residuals) are calculated in any of the
four mentioned ways. Furthermore, we back this up mathematically by showing
that the two correlation measures agree (if the number of SNPs is large) under
the correct admixture model for PCA 1 and PCA 2. For the latter, a few
additional assumptions are required. The estimated covariance (and correlation
coefficient) under the proposed model might be seen as a correction term for
population structure. Subtracting it from the empirical covariance, thus gives
a covariance estimate with baseline zero under the correct model, independent
of the population structure. It is natural to suspect that similar can be done
in models with population structure and kinship, which we will pursue in a
subsequent study.

In the next section, we describe the model, the statistical approach to
compute the residuals, and how we evaluate model fit. In addition, we give
mathematical statements that show how the method performs theoretically. In
the ‘Results’ section, we provide analysis of simulated and real data,
respectively. We end with a discussion. Mathematical proofs are collected in
the appendix.

##  2 Materials and methods

###  2.1 Notation

For an $\ell_{1}\times\ell_{2}$ matrix $A=(A_{ij})_{i,j}$, $A_{\star i}$
denotes the $i$-th column of $A$, $A_{i\star}$ the $i$-th row, $A^{T}$ the
transpose matrix, and $\text{rank}(A)$ the rank. The Frobenius norm of a
square $\ell\times\ell$ matrix $A$ is

$\|A\|_{F}=\sqrt{\sum_{i=1}^{\ell}\sum_{j=1}^{\ell}A_{ij}^{2}}.$

A square matrix $A$ is an orthogonal projection if $A^{2}=A$ and $A^{T}=A$. A
symmetric matrix has $n$ real eigenvalues (with multiplicity) and the
eigenvectors can be chosen such that they are orthogonal to each other. If the
matrix is positive (semi-)definite, then the eigenvalues are positive (non-
negative).

For a random variable/vector/matrix $X$, its expectation is denoted
$\operatorname{\mathbb{E}}[X]$ (provided it exist). The variance of a random
variable $X$ is denoted $\operatorname{var}(X)$, and covariance between two
random variables $X,Y$ is denoted $\operatorname{cov}(X,Y)$ (provided they
exist). Similarly, for a random vector $X=(X_{1},\ldots,X_{n})$, the
covariance matrix is denoted $\operatorname{cov}(X)$. For a sequence $X_{m}$,
$m=0,\ldots,$ of random variables/vectors/matrices, if $X_{m}\to X_{0}$ as
$m\to\infty$ almost surely (convergence for all realisations but a set of zero
probability), we leave out ‘almost surely’ and write $X_{m}\to X_{0}$ as
$m\to\infty$ for convenience.

###  2.2 The PCA and the admixture model

We consider a model with genotype observations from $n$ individuals, and $m$
biallelic sites (SNPs), where $m$ is assumed to be (much) larger than $n$,
$m\geq n$. The genotype $G_{si}$ of SNP $s$ in individual $i$ is assumed to be
a binomial random variable

$G_{si}\sim\operatorname{binomial}(2,\Pi_{si}).$

In matrix notation, we have $G\sim\operatorname{binomial}(2,\Pi)$ with
expectation $\operatorname{\mathbb{E}}(G\mid\Pi)=2\Pi$, where $G$ and $\Pi$
are $m\times n$ dimensional matrices. Conditional on $\Pi$, we assume the
entries of $G$ are independent random variables.

Furthermore, we assume the matrix $\Pi$ takes the form $\Pi=FQ$, where $Q$ is
a (possibly unconstrained) $k\operatorname{\times}n$ matrix of rank $k\leq n$,
and $F$ is a (possibly unconstrained) $m\times k$ matrix, also of rank $k$
(implying $\Pi$ likewise is of rank $k$, Lemma 13). Entry-wise, this amounts
to

$\Pi_{si}=(FQ)_{si}=\sum_{k=1}^{k}F_{sk}Q_{ki},\quad s=1,\dots,n,\quad
i=1,\ldots,m.$

For the binomial assumption to make sense, we must require the entries of
$\Pi$ to be between zero and one.

In the literature, this model is typically encountered in the form of an
admixture model with $k$ ancestral populations, see for example, Pritchard et
al. (2000); Garcia-Erill and Albrechtsen (2020). The general unconstrained
setting which applies to PCA has also been discussed (Cabreros and Storey,
2019). In the case of an admixture model, $Q$ is a matrix of ancestral
admixture proportions, such that the proportion of individual $i$’s genome
originating from population $j$ is $Q_{ji}$. Furthermore, $F$ is a matrix of
ancestral SNP frequencies, such that the frequency of the reference allele of
SNP $s$ in population $j$ is $F_{sj}$. In many applications, the columns of
$Q$ sum to one.

While we lean towards an interpretation in terms of ancestral population
proportions and SNP frequencies, our approach does not enforce or assume the
columns of $Q$ (the admixture proportions) to sum to one, but allow these to
be unconstrained. This is advantageous for at least two reasons. First, a
proposed model might only contain the major ancestral populations, leaving out
older or lesser defined populations. Hence, the sum of ancestral proportions
might be smaller than one. Secondly, when fitting a model with fewer ancestral
populations than the true model, one should only require the admixture
proportions to sum to at most one.

###  2.3 The residuals

Our goal is to design a strategy to assess the hypothesis that $\Pi$ is a
product of two matrices. As we do not know the true $k$, we suggest a number
$k^{\prime}$ of ancestral populations and estimate the model parameters under
this constraint. That is, we assume a model of the form

$G\sim\operatorname{binomial}(2,\Pi_{k^{\prime}}),\quad\Pi_{k^{\prime}}=F_{k^{\prime}}Q_{k^{\prime}},$

where each entry of $G$ follows a binomial distribution. $Q_{k^{\prime}}$ has
dimension $k^{\prime}\operatorname{\times}n$, $F_{k^{\prime}}$ has dimension
$m\operatorname{\times}k^{\prime}$, and
$\text{rank}(Q_{k^{\prime}})=\text{rank}(F_{k^{\prime}})=k^{\prime}$, hence
also $\text{rank}(\Pi_{k^{\prime}})=k^{\prime}$. Throughout, we use the index
$k^{\prime}$ to indicate the imposed rank condition, and assume
$k^{\prime}\leq k$ unless otherwise stated. The latter assumption is only to
guarantee the mathematical validity of certain statements, and is not required
for practical use of the method.

Our approach is build on the residuals, the difference between observed and
predicted data. To define the residuals, we let
$P\colon\operatorname{\mathbb{R}}^{n}\to\operatorname{\mathbb{R}}^{n}$ be the
orthogonal projection onto the $k$-dimensional subspace spanned by the $k$
rows of (the true) $Q$, hence $P=Q^{T}(QQ^{T})^{-1}Q$, and $QP=Q$. Let
$\widehat{P}_{k^{\prime}}$ be an estimate of $P$ based on the data $G$, and
assume $\widehat{P}_{k^{\prime}}$ is an orthogonal projection onto a
$k^{\prime}$-dimensional subspace. Later in this section, we show how an
estimate $\widehat{P}_{k^{\prime}}$ can be obtained from an estimate of
$Q_{k^{\prime}}$ or an estimate of $\Pi_{k^{\prime}}$. Estimates of these
parameters might be obtained using existing methods, based on for example,
maximum likelihood analysis (Wang, 2003; Alexander et al., 2009; Garcia-Erill
and Albrechtsen, 2020). Furthermore, for the three PCA approaches, an estimate
of the projection matrix can simply be obtained from eigenvectors of a
singular value decomposition (SVD) of the data matrix.

We define the $m\times n$ matrix of residuals by

$R_{k^{\prime}}=G-2\widehat{\Pi}=G(I-\widehat{P}_{k^{\prime}}),$

where $G$ is the observed data and $G\widehat{P}_{k^{\prime}}$, the predicted
values. The latter might also be considered an estimate of $2\Pi$, the
expected value of $G$. This definition of residuals is in line with how the
residuals are defined in a multilinear regression model as the difference
between the observed data (here, $G$) and the projection of the data onto the
subspace spanned by the regressors (here, $G\widehat{P}_{k^{\prime}}$). The
essential difference being that in a multilinear regression model, the
regressors are known and does not depend on the observed data, while
$\widehat{P}_{k^{\prime}}$ is estimated from the data.

We assess the model fit by studying the correlation matrix of the residuals in
two ways. First, we consider the _empirical covariance matrix_ $\widehat{B}$
with entries

$\displaystyle\widehat{B}_{ij}$
$\displaystyle=\frac{1}{m-1}\sum_{s=1}^{m}(R_{k^{\prime},si}-{\hbox{\vbox{\hrule
height=0.5pt\kern 1.29167pt\hbox{\kern-1.00006pt$R$\kern
0.0pt}}}}_{k^{\prime},i})(R_{k^{\prime},sj}-{\hbox{\vbox{\hrule
height=0.5pt\kern 1.29167pt\hbox{\kern-1.00006pt$R$\kern
0.0pt}}}}_{k^{\prime},j})$
$\displaystyle=\frac{1}{m-1}\sum_{s=1}^{m}(R_{k^{\prime},si}R_{k^{\prime},sj}-{\hbox{\vbox{\hrule
height=0.5pt\kern 1.29167pt\hbox{\kern-1.00006pt$R$\kern
0.0pt}}}}_{k^{\prime},i}\ {\hbox{\vbox{\hrule height=0.5pt\kern
1.29167pt\hbox{\kern-1.00006pt$R$\kern 0.0pt}}}}_{k^{\prime},j}),$

where

${\hbox{\vbox{\hrule height=0.5pt\kern 1.29167pt\hbox{\kern-1.00006pt$R$\kern
0.0pt}}}}_{k^{\prime}i}=\frac{1}{m}\sum_{s=1}^{m}R_{k^{\prime},si},$

and the corresponding _empirical correlation matrix_ with entries

$\widehat{b}_{ij}=\frac{\widehat{B}_{ij}}{\sqrt{\widehat{B}_{ii}\widehat{B}_{jj}}},$

$i,j=1,\ldots,n$. Secondly, we consider the _estimated covariance matrix_

$\widehat{C}=(I-\widehat{P}_{k^{\prime}})\widehat{D}(I-\widehat{P}_{k^{\prime}})$

with corresponding _estimated correlation matrix_ ,

$\displaystyle\widehat{c}_{ij}=\frac{\widehat{C}_{ij}}{\sqrt{\widehat{C}_{ii}\widehat{C}_{jj}}},$

$i,j=1,\ldots,n$. Here, $\widehat{D}$ is the $n\times n$ diagonal matrix
containing the average heterozygosities of each individual,

$\widehat{D}_{ii}=\frac{1}{m}\sum_{s=1}^{m}G_{si}(2-G_{si}),\quad
i=1,\ldots,n.$

Under reasonable regularity conditions, we can quantify the behaviour of
$\widehat{B}$ and $\widehat{C}$ as the number of SNPs become large.
Specifically, we assume the rows of $F$ are independent and identically
distributed with distribution $\text{Dist}(\mu,\Sigma)$, where $\mu$ denote
the $k$-dimensional mean vector of the distribution, and $\Sigma$ the
$k\operatorname{\times}k$-covariance matrix, that is,

$\displaystyle F_{s\star}=(F_{s1},\ldots,F_{sk})\,$
$\displaystyle\stackrel{{\scriptstyle\text{iid}}}{{\sim}}\,\text{Dist}(\mu,\Sigma),$

$s=1,\ldots,m$. The matrix $Q$ is assumed to be non-random, that is, fixed.
These assumptions are standard and typically used in simulation of genetic
data, see for example, Pickrell and Pritchard (2012); Cabreros and Storey
(2019); Garcia-Erill and Albrechtsen (2020). Often $\text{dist}(\mu,\Sigma)$
is taken to be the product of $k$ independent uniform distributions in which
case $\mu=0.5(1,1,\ldots,1)$ and $\Sigma$ is a diagonal matrix with entries
$1/12$, though other choices have been applied, see for example Balding and
Nichols (1995); Conomos et al. (2016).

Let $D$ be the diagonal matrix with entries

$D_{ii}=2\operatorname{\mathbb{E}}[\Pi_{si}(1-\Pi_{si})],\quad i=1,\ldots,n.$
(1)

It follows from Lemma 7 in the appendix, that $\widehat{D}$ converges to $D$
as $m\to\infty$. Furthermore, as $D_{ii}$ is the variance of $G_{si}$ (it is
binomial), then $\widehat{D}_{ii}$ might be considered an estimate of this
variance. The proofs of the statements are in the appendix.

######  Theorem 1.

Let $k^{\prime}\leq k$. Under the given assumptions, suppose further that
$\widehat{P}_{k^{\prime}}\to P_{k^{\prime}}$ as $m\to\infty$, for some matrix
$P_{k^{\prime}}$. Then, $P_{k^{\prime}}$ is an orthogonal projection.
Furthermore, the following holds,

$\displaystyle\widehat{B}$
$\displaystyle\,\to\,(I-P_{k^{\prime}})(D+4Q^{T}\Sigma Q)(I-P_{k^{\prime}}),$
$\displaystyle\widehat{C}$
$\displaystyle\,\to\,(I-P_{k^{\prime}})D(I-P_{k^{\prime}}),$

as $m\to\infty$. Hence, also

$\displaystyle\widehat{B}\,-\,\widehat{C}\,\,$
$\displaystyle\to\,\,4(I-P_{k^{\prime}})Q^{T}\Sigma Q(I-P_{k^{\prime}})$
$\displaystyle\,\,=\,\,4(P-P_{k^{\prime}})Q^{T}\Sigma Q(P-P_{k^{\prime}}),$

as $m\to\infty$. For $k^{\prime}=k$, if $P_{k}=P$, then the right hand side is
the zero matrix, whereas this is not the case in general for $k^{\prime}<k$.

######  Theorem 2.

Assume $k^{\prime}=k$ and $P_{k}=P$. Furthermore, suppose as in Theorem 1 and
that the vector with all entries equal to one is in the space spanned by the
rows of $Q$ (this is, for example, the case if the admixture proportions sum
to one for each individual). Then,

$\displaystyle\frac{\sum_{i=1}^{n}\sum_{j=1,i\not=j}^{n}\widehat{B}_{ij}}{\sum_{i=1}^{n}\widehat{B}_{ii}}$
$\displaystyle\to\ -1,\quad\text{as}\quad m\to\infty.$ (2)

In addition, if $Q$ takes the form

$Q=\begin{pmatrix}Q_{1}&0&\cdots&0\\\ 0&Q_{2}&\cdots&0\\\
\vdots&\vdots&\ddots&\vdots\\\ 0&0&\cdots&Q_{r}\end{pmatrix}$

where $Q_{\ell}$ has dimension $k_{\ell}\times n_{\ell}$,
$\sum_{\ell=1}^{r}k_{\ell}=k$ and $\sum_{\ell=1}^{r}n_{\ell}=n$, then (2)
holds for each component of $n_{\ell}$ individuals. If $Q_{\ell}=(1\ldots 1)$,
then

$\widehat{b}_{ij}\ \to\ -\frac{1}{n_{\ell}-1},\quad\text{as}\quad m\to\infty,$

for all individuals $i,j$ in the $\ell$-th component, irrespective the form of
$Q_{\ell^{\prime}}$, $\ell^{\prime}\not=\ell$.

######  Theorem 3.

Assume $k^{\prime}=k$ and $P_{k}=P$. Furthermore, suppose as in Theorem 1 and
that $Q$ takes the form

$Q=\begin{pmatrix}Q_{1}&Q_{2}\\\ 0&Q_{3}\end{pmatrix},$

where $Q_{1}=(1\ldots 1)$ has dimension $1\times n_{1}$, $n_{1}\leq n$. Then,
$\widehat{b}_{ij}$ converges as $m\to\infty$ to a value larger than or equal
to $-\tfrac{1}{n_{1}-1},$ for all $i,j=1,\ldots,n_{1}$.

The same statements in the last two theorems hold with $\widehat{B}$ and
$\widehat{b}$ replaced by $\widehat{C}$ and $\widehat{c}$, respectively.

The three theorems provide means to evaluate the model. In particular, Theorem
1 might be used to assess the correctness (or appropriateness) of the proposed
$k^{\prime}$, while Theorem 2 and Theorem 3 might be used to assess whether
data from a group of individuals (e.g., a modern day population) originates
from a single ancestral population, irrespective, the origin of the remaining
individuals. We give examples in the Results section.

The work flow is shown in Algorithm 1. We process real and simulated genotype
data using PCA 1, PCA 2, PCA 3, and the software ADMIXTURE, and evaluate the
fit of the model.

  1. 1.

Choose $k^{\prime}$,

  2. 2.

Compute an estimate $\widehat{P}_{k^{\prime}}$ of the projection $P$,

  3. 3.

Calculate the residuals $R_{k^{\prime}}=G(I-\widehat{P}_{k^{\prime}})$,

  4. 4.

Calculate the correlation coefficients, $\widehat{b}$ and $\widehat{c}$,

  5. 5.

Plot $\widehat{b}$ and the difference, the corrected correlation coefficients,
$\widehat{b}-\widehat{c}$,

  6. 6.

Assess visually the fit of the model.

Algorithm 1 Work flow of the proposed method

###  2.4 Estimation of $P_{k^{\prime}}$

Estimation of $Q,F,$ and $\Pi$ has received considerable interest in the
literature, using for example, maximum likelihood (Wang, 2003; Alexander et
al., 2009), Bayesian approaches (Pritchard et al., 2000) or PCA (Engelhardt
and Stephens, 2010).

We discuss different ways to obtain an estimate $\widehat{P}_{k^{\prime}}$ of
$P$.

####  2.4.1 Using an estimate $\widehat{Q}_{k^{\prime}}$ of $Q_{k^{\prime}}$

An estimate $\widehat{P}_{k^{\prime}}$ might be obtained by projecting onto
the subspace spanned by the $k^{\prime}$ rows of $\widehat{Q}_{k^{\prime}}$,

$\widehat{P}_{k^{\prime}}=\widehat{Q}_{k^{\prime}}^{T}(\widehat{Q}_{k^{\prime}}\widehat{Q}_{k^{\prime}}^{T})^{-1}\widehat{Q}_{k^{\prime}},$

assuming $\text{rank}(\widehat{Q}_{k^{\prime}})=k^{\prime}$ for the
calculation to be valid.

We apply this approach to estimate the projection matrix using output from the
software ADMIXTURE.

####  2.4.2 Using an estimate $\widehat{\Pi}_{k^{\prime}}$ of
$\Pi_{k^{\prime}}$

Let $\widetilde{\Pi}_{k^{\prime}}$ be $k^{\prime}$ linearly independent rows
chosen from $\widehat{\Pi}_{k^{\prime}}$ (out of $m$ rows). Then, an estimate
$\widehat{P}_{k^{\prime}}$ of $P_{k^{\prime}}$ is

$\widehat{P}_{k^{\prime}}=\widetilde{\Pi}_{k^{\prime}}^{T}(\widetilde{\Pi}_{k^{\prime}}\widetilde{\Pi}_{k^{\prime}}^{T})^{-1}\widetilde{\Pi}_{k^{\prime}},$

assuming $\text{rank}(\widehat{\Pi}_{k^{\prime}})=k^{\prime}$ for the
calculation to be valid. Alternatively, one might apply the Gram-Schmidt
method in which case the vectors are orthonormal by construction and
$\widehat{P}_{k^{\prime}}=\widetilde{\Pi}_{k^{\prime}}^{T}\widetilde{\Pi}_{k^{\prime}}$.
The estimate $\widehat{P}_{k^{\prime}}$ is independent of the choice of the
$k^{\prime}$ rows, provided
$\text{rank}(\widehat{\Pi}_{k^{\prime}})=k^{\prime}$.

####  2.4.3 Using PCA 1

We consider a PCA approach, originally due to Chen and Storey (2015), to
estimate the space spanned by the rows of $Q$. We follow the procedure laid
out in Cabreros and Storey (2019).

Let $\widehat{H}$ be the symmetric matrix

$\widehat{H}=\frac{1}{m}G^{T}G-\widehat{D}.$

Since $\widehat{H}$ is symmetric, all eigenvalues are real and the matrix is
diagonalisable. Furthermore, $\widehat{H}$ is a variance adjusted version of
$\frac{1}{m}G^{T}G$, see (1). Let $u_{1},\ldots,u_{k^{\prime}}$ be
$k^{\prime}\leq k$ orthogonal eigenvectors belonging to the $k^{\prime}$
largest eigenvalues of $\widehat{H}$, counted with multiplicities. Define the
$n\operatorname{\times}k^{\prime}$ matrix
$U_{k^{\prime}}=(u_{1},\ldots,u_{k^{\prime}})$ and the
$n\operatorname{\times}n$ orthogonal projection matrix

$\widehat{P}_{k^{\prime}}=U_{k^{\prime}}(U_{k^{\prime}}^{T}U_{k^{\prime}})^{-1}U_{k^{\prime}}^{T}=U_{k^{\prime}}U_{k^{\prime}}^{T}$

onto the subspace given by the span of the vectors
$u_{1},\ldots,u_{k^{\prime}}$.

In this particular case, convergence of $\widehat{P}_{k^{\prime}}$ can be made
precise. Define the matrix $H=4Q^{T}(\Sigma+\mu\mu^{T})Q$. Then, $H$ is
symmetric and positive semi-definite because $\Sigma$ and $\mu\mu^{T}$ both
are positive semi-definite. Hence, $H$ has non-negative eigenvalues.
Furthermore, according to Lemma 8 in the appendix, $\widehat{H}$ converges to
$H$ as $m\to\infty$.

######  Theorem 4.

Assume $k^{\prime}\leq k$. Let $\lambda_{1}\geq\ldots\geq\lambda_{n}\geq 0$ be
the eigenvalues of $H$, with corresponding orthogonal eigenvectors
$v_{1},\ldots,v_{n}$. In particular, $\lambda_{k+1}=\ldots=\lambda_{n}=0$, as
$Q$ has rank $k$. Let $P_{k^{\prime}}$ be the orthogonal projection onto the
span of $v_{1},\ldots,v_{k^{\prime}}$, that is,

$P_{k^{\prime}}=V_{k^{\prime}}(V_{k^{\prime}}^{T}V_{k^{\prime}})^{-1}V_{k^{\prime}}^{T}=V_{k^{\prime}}V_{k^{\prime}}^{T},$

where $V_{k^{\prime}}=(v_{1},\ldots,v_{k^{\prime}})$.

Assume $k^{\prime}=n$ or $\lambda_{k^{\prime}}>\lambda_{k^{\prime}+1}$,
referred to as the eigenvalue condition. Then, $\widehat{P}_{k^{\prime}}\to
P_{k^{\prime}}$ as $m\to\infty$. If the eigenvalue condition is fulfilled for
$k^{\prime}=k$, then $P_{k}=P$, that is, $P_{k}$ is the orthogonal projection
onto the span of the row vectors of $Q$. In particular, the eigenvalue
condition is fulfilled for $k^{\prime}=k$ if and only if $\Sigma+\mu\mu^{T}$
is positive definite. The latter is the case if $\Sigma$ is positive definite.

For $k^{\prime}=k$, the correct row space of $Q$ is found eventually, but not
$Q$ itself. If $k^{\prime}<k$, then a subspace of this row space is found,
corresponding to the $k^{\prime}$ largest eigenvalues. As the data is not mean
centred, we discard the first principal component, and use the subsequent
$k^{\prime}-1$ eigenvectors and eigenvalues.

####  2.4.4 Using PCA 2 (mean centred data)

A popular approach to estimation of $\Pi$ in the admixture model is PCA based
on mean centred data, or mean and variance normalised data (Pritchard et al.,
2000; Engelhardt and Stephens, 2010; Patterson et al., 2006).

Let $G_{1}=G-\tfrac{1}{n}GE=G(I-\tfrac{1}{n}E)$ be the SNP-wise mean centred
genotypes, where $E$ is an $n\times n$ matrix with all entries equal to one.
Following the exposition and notation in Cabreros and Storey (2019), let
$G_{1}=U\Delta V^{T}$ be the SVD of $G_{1}$, where $\Delta V^{T}$ consists of
the row-wise principal components of $G_{1}$, ordered according to the
singular values. Define

$S_{k^{\prime}}=\begin{pmatrix}U^{T}_{1:(k^{\prime}-1)}\\\ e\end{pmatrix},$

where $e=(1\,1\,\ldots\,1)$ is a vector with all entries one, and
$U^{T}_{1:(k^{\prime}-1)}$ contains the top $k^{\prime}-1$ rows of $U^{T}$.
Then, an estimate of the projection is

$\widehat{P}_{k^{\prime}}=S_{k^{\prime}}^{T}(S_{k^{\prime}}S_{k^{\prime}}^{T})^{-1}S_{k^{\prime}}.$

The squared singular values in the SVD decomposition of $G_{1}$ are the same
as the eigenvalues of

$\widehat{H}_{1}=\frac{1}{m}G^{T}_{1}G_{1}=\frac{1}{m}\left(I-\frac{1}{n}E\right)G^{T}G\left(I-\frac{1}{n}E\right)$

(Jolliffe, 2002). We have

$\displaystyle\operatorname{\mathbb{E}}[\widehat{H}_{1}]$
$\displaystyle=\frac{1}{m}\left(I-\frac{1}{n}E\right)\operatorname{\mathbb{E}}[G^{T}G]\left(I-\frac{1}{n}E\right)$
$\displaystyle=\left(I-\frac{1}{n}E\right)(D+4Q^{T}(\Sigma+\mu\mu^{T})Q)\left(I-\frac{1}{n}E\right).$
(3)

Let $H_{1}$ denote the right hand side of (3).

######  Theorem 5.

Let $\lambda_{1}\geq\ldots\geq\lambda_{n}$ be the eigenvalues of $H_{1}$, with
corresponding orthogonal eigenvectors $v_{1},\ldots,v_{n}$. In particular,
$v_{n}=e$ and $\lambda_{n}=0$. If $D$ has all diagonal entries positive, then
$\lambda_{n-1}>0$.

Let $k^{\prime}\leq n$ and let $P_{k^{\prime}}$ be the orthogonal projection
onto the span of $v_{1},\ldots,v_{k^{\prime}-1},e$, that is,

$P_{k^{\prime}}=V_{k^{\prime}}(V_{k^{\prime}}^{T}V_{k^{\prime}})^{-1}V_{k^{\prime}}^{T},$

where $V_{k^{\prime}}=(v_{1},\ldots,v_{k^{\prime}-1},e)$. If $k^{\prime}=n$ or
$\lambda_{k^{\prime}}>\lambda_{k^{\prime}+1}$, then
$\widehat{P}_{k^{\prime}}\to P_{k^{\prime}}$ as $m\to\infty$.

There are no guarantees that for $k^{\prime}=k$, we have $P_{k}=P$ and that
the difference between $\widehat{B}$ and $\widehat{C}$ converges to zero for
large $m$. However, this is the case under some extra conditions, and appears
to be the case in many practical situations, see the Results section.

######  Theorem 6.

Assume $D=dI$ for some $d>0$. Furthermore, assume the vector $e$ is in the row
space of $Q$ (this is, for example, the case if the admixture proportions sum
to one for each individual). Then, $\lambda_{k}=\ldots=\lambda_{n-1}=d$, and
$\lambda_{n}=0$.

If $\Sigma+\mu\mu^{T}$ is positive definite, then $\lambda_{k+1}>\lambda_{k}$
and $P_{k}=P$, where $P_{k}$ is as in Theorem 4. As a consequence, with
$k^{\prime}=k$ in Theorem 1, $\widehat{B}-\widehat{C}\to 0$ as $m\to\infty$.

####  2.4.5 Using PCA 3 (mean and variance normalised data)

Let $G_{2}=W^{-1}G_{1}$ be the SNP mean and variance normalised genotypes,
where $W$ is an $m^{\prime}\times m^{\prime}$ diagonal matrix with $s$-th
entry being the observed standard deviation of the genotypes of SNP $s$. All
SNPs for which no variation are observed are removed, hence the number of SNPs
might be smaller than the original number, $m^{\prime}\leq m$. Following the
same procedure as for PCA 2, let $G_{2}=U\Delta V^{T}$ be the SVD of $G_{2}$,
where $\Delta V^{T}$ consists of the row-wise principal components of $G_{2}$,
ordered according to the singular values. Define

$S_{k^{\prime}}=\begin{pmatrix}V^{T}_{1:(k^{\prime}-1)}\\\ e\end{pmatrix},$

where $e=(1\,1\,\ldots\,1)$, and $V^{T}_{1:(k^{\prime}-1)}$ contains the top
$k^{\prime}-1$ rows of $V^{T}$. Then, an estimate of the projection is
$\widehat{P}_{k^{\prime}}=S_{k^{\prime}}^{T}(S_{k^{\prime}}S_{k^{\prime}}^{T})^{-1}S_{k^{\prime}}.$

We are not aware of any theoretical justification of this procedure similar to
Theorem 1, but it appears to perform well in many practical situations,
according to our simulations.

###  2.5 Simulation of genotype data

We simulated genotype data from different demographic scenarios using
different sampling strategies. We deliberately choose different sampling
strategies to challenge the method. We first made simple simulations that
illustrate the problem of model fit as well as to demonstrate the theoretical
and practical properties of the residual correlations that arise from having
data from a finite number of individuals and a large number of SNPs. An
overview of the simulations are given in Table 1.

In the first two scenarios, the ancestral allele frequencies are simulated
independently for each ancestral population from a uniform distribution,
$F_{si}\sim\text{Unif}(0,1)$ for each site $s=1,\ldots,m$ and each ancestral
population $i=1,\ldots,k$. In scenario 1, we simulated unadmixed individuals
from three populations with either an equal or an unequal number of sampled
individuals from each population. In scenario 2, we simulated two ancestral
populations and a population that is admixed with half of its ancestry coming
from each of the two ancestral populations.

In scenario 3, we set $F_{si}\sim\text{Unif}(0.01,0.99)$ and simulated spatial
admixture in a way that resembles a spatial decline of continuous gene flow
between populations living in a long narrow island. We first simulated a
single population in the middle of the long island. From both sides of the
island, we then recursively simulated new populations from a Balding-Nichols
distribution with parameter $F_{st}=0.001$ using the R package ‘bnpsd’ (Ochoa
and Storey, 2019). In this way, each pair of adjacent populations along the
island has an $F_{st}$ of 0.001. Additional details on the simulation and an
schematic visualization can be found in Figure 2 of Garcia-Erill and
Albrechtsen (2020).

In scenario 4, we first simulated allele frequencies for an ancestral
population from a symmetric beta distribution with shape parameter 0.03,
$F_{si}\sim\text{Beta}(0.3,0.3)$, which results in an allele frequency
spectrum enriched for rare variants, mimicking the human allele frequency
spectrum. We then sampled allele frequencies from a bifurcating tree
(((pop1:0.1,popGhost:0.2):0.05,pop2:0.3):0.1,pop3:0.5), where pop1 and
popGhost are sister populations and pop3 is an outgroup. Using the Balding-
Nichols distribution and the $F_{st}$ branch lengths of the tree (see Figure
5), we sampled allele frequencies in the four leaf nodes. Then, we created an
admixed population with 30% ancestry from popGhost and 70% from pop2. We
sampled 10 million genotypes for 50 individuals from each population except
for the ghost population which was not included in the analysis, and
subsequently removed sites with a sample minor allele frequency below 0.05,
resulting in a total of 694,285 sites.

In scenario 5, we simulated an ancestral population with allele frequencies
from a uniform distribution $F_{si}\sim$ Unif $(0.05,0.95)$, from which we
sampled allele frequencies for two daughter populations from a Balding
Nicholds distributions with $F_{st}=0.3$ from the ancestral population, using
’bnpsd’. We then created recent hybrids based on a pedigree where all but one
founder has ancestry from the first population. The number of generations in
the pedigree then determines the admixture proportions and the age of the
admixture where F1 individuals have one unadmixed parent from each population
and backcross individuals have one unadmixed parent and the other F1. Double
backcross individuals have one unadmixed parent and the other is a backcross.
We continue to quadruple backcross with one unadmixed parent and the other
triple backcross. Note that for the recent hybrids the ancestry of the pair of
alleles at each loci is no longer independent which is a violation of the
admixture model.

Table 1: Overview of simulations.

Scenario | $\bm{k}$ | $\bm{n}$ | $\bm{m}$ | Description |  $\bm{F}_{is}$111Ancestral allele frequencies, $i=1,\ldots,k$
---|---|---|---|---|---
1 | 3 | 20,20,20 | $500K$ | Unadmixed | $\text{Unif}(0,1)$
1 | 3 | 10,20,30 | $500K$ | Unadmixed | $\text{Unif}(0,1)$
2 | 2 | 20,20,20 | $500K$ | Admixed | $\text{Unif}(0,1)$
2 | 2 | 10,20,30 | $500K$ | Admixed | $\text{Unif}(0,1)$
3 |  | 500 |  $100K$222after applying MAF$>5$% filtering, 88,082 remained. | Spatial with $F_{st}=0.001$ | $\text{Unif}(0.01,0.99)$
|  |  |  | between adjacent populations |
4 | 4 | 50,50,50,50,0333No reference samples are provided on the ghost population. |  $10M$444after applying MAF$>5$% filtering, 694,285 remained. | Ghost admixture | $\text{Beta}(0.3,0.3)$
5 | 2 | 20,20,50 | $500K$ | Recent hybrids | $\text{Unif}(0.05,0.95)$

##  3 Results

###  3.1 Scenario 1

In this first set-up, we demonstrate the method using PCA 1 only. We simulated
unadmixed individuals from $k=3$ ancestral populations

$Q=\begin{pmatrix}{1}_{n_{1}}&0&0\\\ 0&{1}_{n_{2}}&0\\\
0&0&{1}_{n_{3}}\end{pmatrix},$

where ${1}_{n_{i}}$ is a row vector with all elements being one, and
$n_{1}+n_{2}+n_{3}=n$. We simulated genotypes for $n=60$ individuals with
sample sizes $n_{1},n_{2}$ and $n_{3}$, respectively, as detailed in the
previous section. In Figure 1(A), we show the residual correlation
coefficients for $k^{\prime}=2,3$ and plot the corresponding major PCs. For
the PCA 1 approach, the first principal component does not relate to
population structure as the data is not mean centered, and we use the
following $k^{\prime}-1$ principal components.

When assuming that there are only two populations, $k^{\prime}=2$, we note
that the empirical correlation coefficients appear largely consistent within
each population sample, but the corrected correlation coefficients are
generally non-zero with different signs, which points to model misfit. In
contrast, when assuming the correct number of populations is $k^{\prime}=3$,
the empirical correlation coefficients match nicely the theoretical values of
$-\tfrac{1}{n_{i}-1}$, which comply with Theorem 2 (see Table 2). A fairly
homogeneous pattern in the corrected correlation coefficients appears around
zero across all samples. This is a good indication that the model fits well
and that the PCA plots using principal components 2 and 3 reflex the data
well.

Figure 1: Results for simulated Scenario 1. (A) The upper triangle in the
plots shows the empirical correlation coefficients $\hat{b}$ and the lower
triangle shows the corrected correlation coefficients $\hat{b}-\hat{c}$. (B)
The major principal components ($k^{\prime}=3$) result in a clear separation
of the three samples (all data points within each sample are almost
identical).  Table 2: The mean (standard deviation) of $\hat{b}$ and
$\hat{b}-\hat{c}$ within each population using PCA 1.

Scenario 1 | $\bm{k}^{\prime}$ | $\bm{n}$ |  | pop1 | pop2 | pop3 |
---|---|---|---|---|---|---|---
| $3$ | $(20,20,20)$ |  $\hat{b}\ $555The second line of $\hat{b}$ in each case shows the theoretical value obtained from the limit in Theorem 1. | -0.0526 (0.0015) | -0.0526 (0.0016) | -0.0526 (0.0016) |
|  |  |  | -0.0526 | -0.0526 | -0.0526 |
|  |  | $\hat{b}-\hat{c}$ | 0e-04 (0.0015) | 0e-04 (0.0016) | 0e-04 (0.0016) |
|  | $(10,20,30)$ | $\hat{b}$ | -0.1111 (0.0011) | -0.0526 (0.0016) | -0.0345 (0.0016) |
|  |  |  | -0.1111 | -0.0526 | -0.0345 |
|  |  | $\hat{b}-\hat{c}$ | 0e-04 (0.0012) | 0e-04 (0.0016) | 0e-04 (0.0016) |
Scenario 2 | $\bm{k}^{\prime}$ | $\bm{n}$ |  | pop1 | admixed | pop3 |
| $2$ | $(20,20,20)$ | $\hat{b}$ | -0.0419 (0.0015) | -0.0192 (0.0015) | -0.0420 (0.0015) |
|  |  |  | -0.0420 | -0.0193 | -0.0420 |
|  |  | $\hat{b}-\hat{c}$ | 0e-04 (0.0015) | 0e-04 (0.0015) | 0e-04 (0.0015) |
|  | $(10,20,30)$ | $\hat{b}$ | -0.0701 (0.0018) | -0.0228 (0.0014) | -0.0304 (0.0016) |
|  |  |  | -0.0701 | -0.0229 | -0.0304 |
|  |  | $\hat{b}-\hat{c}$ | 0e-04 (0.0017) | 0e-04 (0.0014) | 0e-04 (0.0016) |
Scenario 4 | $\bm{k^{\prime}}$ | $\bm{n}$ |  | pop1 | pop2 | pop3 | pop4
| $3$ | $(50,50,50,50)$ | $\hat{b}$ | -0.0190 (0.0015) | 0.0027 (0.0015) | -0.0204 (0.0017) | 0.0122 (0.0013)
|  |  | $\hat{b}-\hat{c}$ | 0.0009 (0.0015) | 0.0147 (0.0015) | 0e-04 (0.0017) | 0.0208 (0.0013)
| $4$ |  | $\hat{b}$ | -0.0204 (0.0015) | -0.0204 (0.0015) | -0.0204 (0.0017) | -0.0204 (0.0014)
|  |  | $\hat{b}-\hat{c}$ | 0e-04 (0.0015) | 0e-04 (0.0015) | 0e-04 (0.0017) | 0e-04 (0.0013)

###  3.2 Scenario 2

In this set-up we also include admixed individuals. We simulated samples from
two ancestral populations and individuals that are a mix of the two. We then
applied all three PCA procedures and the software ADMIXTURE to the data.
Specifically, we choose

$Q=\begin{pmatrix}{1}_{n_{1}}&\frac{1}{2}{1}_{n_{2}}&0\\\
0&\frac{1}{2}{1}_{n_{2}}&{1}_{n_{3}}\end{pmatrix},$

with $k=2$ true ancestral populations, and $(n_{1},n_{2},n_{3})=(20,20,20)$ or
$(n_{1},n_{2},n_{3})=(10,20,30)$, see the previous section for details. We
analysed the data with $k^{\prime}=1,2,3$, and obtained the correlation
structure shown in Figures 2 and 3, and Table 2. The two standard approaches
PCA 2 and PCA 3 show almost identical results, hence only PCA 2 is shown in
the figures. Both PCA 2 and PCA 3 use the top principal components, while PCA
1 disregards the first, hence the discrepancy in the axis labeling in Figures
2(b) and 3/b). For $k^{\prime}=1$ none of the principal components are used
and the predicted normalized genotypes is simply 0. All four methods show
consistent results, in particular, for the correct $k^{\prime}$ ($=2$), while
there are smaller discrepancies between the methods for wrong
$k^{\prime}=1,3$. This is most pronounced for PCA 1 and ADMIXTURE. We note
that the average correlation coefficient of $\widehat{b}$ within each
population sample comply with Theorem 1 (see Table 2). A fairly homogeneous
pattern in the corrected correlation coefficients appears around zero across
all samples for $k^{\prime}=2$, as in scenario 1, which shows that the model
fits well. However, unlike in scenario 1 the bias for the empirical
correlation coefficient is not a simple function of the sample size (see Table
2).

In this case, and similarly in all other investigated cases, we don’t find any
big discrepancies between the four methods. Therefore, we only show the
results of PCA 1 for which we have theoretical justification for the results.

Figure 2: Results for simulated Scenario 2 with equal sample sizes. (A) For
each of PCA 1, PCA 2 and ADMIXTURE, the upper left triangle in the plots shows
the empirical correlation $\hat{b}$ and the lower right triangle shows the
difference $\hat{b}-\hat{c}$ with sample sizes
$(n_{1},n_{2},n_{3})=(20,20,20)$. (B) The major principal component for the
PCA based methods for $k^{\prime}=2$ (in which case there is only one
principal component). Individuals within each sample have the same color. (C)
The estimated admixture proportions in the case of ADMIXTURE.  Figure 3:
Results for simulated Scenario 2 with unequal sample sizes. (A) For each of
PCA 1, PCA 2 and ADMIXTURE, the upper left triangle in the plots shows the
empirical correlation $\hat{b}$ and the lower right triangle shows the
difference $\hat{b}-\hat{c}$ with sample sizes
$(n_{1},n_{2},n_{3})=(20,20,20)$. (B) The major principal component for the
PCA based methods for $k^{\prime}=2$ (in which case there is only one
principal component). Individuals within each sample have the same color. (C)
The estimated admixture proportions in the case of ADMIXTURE.

###  3.3 Scenario 3

We simulated genotypes for $n=500$ individuals at $m=88,082$ sites with
continuous genetic flow between individuals, thus there is not a true $k$. We
analysed the data assuming $k^{\prime}=2,3$, see Figure 4. In the figure, the
individuals are ordered according to the estimated proportions of the
ancestral populations, hence it appears there is a color wave pattern in the
empirical and the corrected correlation coefficients, see Figure 4(A). As
expected, the corrected correlation coefficients are closer to zero for
$k^{\prime}=3$ than $k^{\prime}=2$, though the deviations from zero are still
large. We thus find no support for the model for either value of $k^{\prime}$.
This is consistent with the plots of the major PCs, that show continuous
change without grouping the data into two or three clusters, see Figure 4(B).

Figure 4: Results for simulated scenario 3. (A) The upper triangle in the
plots shows the empirical correlation $\hat{b}$ and the lower triangle shows
the difference $\hat{b}-\hat{c}$. (B) The major principal components (only one
in the case of $k^{\prime}=2$).

###  3.4 Scenario 4

This case is based on the tree in Figure 5, which include an unsampled (so-
called) ghost population, popGhost. The popGhost is sister population to pop1.

Figure 5: Schematic of the tree used to simulate population allele frequencies
for Scenario 4, including 5 populations: pop1, pop2, pop3, pop4 and popGhost.
The pop4 population is the result of admixture between pop2 and popGhost, for
which there are no individuals sampled and is therefore a ghost population.
The values in the branches indicate the drift in units of $F_{ST}$. The values
along the two admixture edges are the admixture proportions coming from each
population.

We simulated genotypes for $n=200$ individuals: 150 unadmixed samples from
pop1, pop2, and pop3; and 50 samples admixed with 0.3 ancestry from popGhost
and 0.7 ancestry from pop2 (as pop4), as detailed in the previous section. As
there is drift between the populations and hence genetic differences, the
correct $k=4$ (pop1, pop2, pop3, popGhost). This is picked up by our method
that clearly shows $k^{\prime}=3$ is wrong with large deviation from zero in
the corrected correlation coefficients. In contrast, for $k^{\prime}=4$, the
corrected correlation coefficients are almost zero (Figure 6).

Figure 6: Results for simulated scenario 4. (A) The upper triangle in the
plots shows the empirical correlation $\hat{b}$ and the lower triangle shows
the difference $\hat{b}-\hat{c}$. (B) The major principal components for
$k^{\prime}=4$, that result in a clear separation of the four samples (all
data points within each sample are almost identical).

###  3.5 Scenario 5

In the last example, we simulated two populations (originating from a common
ancestral population) and created admixed populations by backcrossing, as
detailed in the previous section. Thus, the model does not fulfil the
assumptions of the admixture model in that the number of reference alleles are
not binomially distributed, but depends on the particular backcross and the
frequencies of the parental populations.

We simulate genotypes for $n=90$ individuals at $m=500,000$ sites. There are
20 homogeneous individuals from each parental population, and 10 different
individuals from each of the different recent admixture classes. Then, we
analysed the data with $k^{\prime}=2$ and found the corrected correlation
coefficients deviated consistently from zero, in particular for one of the
parental populations (Figure 7). We are thus able to say the admixture model
does not provide a reasonable fit.

Figure 7: Results for simulated scenario 5 (recent admixture). (A) The upper
triangle in the plots shows the empirical correlation $\hat{b}$ and the lower
triangle shows the difference $\hat{b}-\hat{c}$. (B) The major principal
component for $k^{\prime}=2$.

###  3.6 Real data

We analysed a whole genome sequencing data set from the 1000 Genomes Project
(Auton et al., 2015), see also Garcia-Erill and Albrechtsen (2020) where the
same data is used. It consists of data from five groups of different descent:
a Yoruba group from Ibadan, Nigeria (YRI), residents from Southwest US with
African ancestry (ASW), Utah residents with Northern and Western European
ancestry (CEU), a group with Mexican ancestry from Los Angeles, California
(MXL), and a group of Han Chinese from Beijing, China (CHB) with sample sizes
$108,61,99,63$ and $103$, respectively, in total, $n=434$. We kept only sites
present in the Human Origins SNP panel (Lazaridis et al., 2014), with a total
of $m=406,279$ SNPs were left after a MAF filter of 0.05.

We analyzed the data with $k^{\prime}=3,4$. For $k^{\prime}=3$, Figure 8 shows
that it is not possible to explain the relationship between MXL, CEU and CHB,
indicating that MXL is not well explained as a mixture of the two. For
$k^{\prime}=4$, the color shades of the corrected correlation coefficients are
almost negligible within each population, pointing at a contribution from a
native american population. This is further corroborated in Figure 8(D) that
shows estimated proportions from the four ancestral populations using the
software ADMIXTURE.

Figure 8: The residual correlation coefficient, the inferred population
structure and the admixture proportions of a real human data from 1000 Genomes
project. (A) The upper triangle in the plots shows the empirical correlation
coefficient $\hat{b}$ and the lower triangle shows the difference
$\hat{b}-\hat{c}$. (B) The three major principal component for PCA 1 for
$k^{\prime}=4$. (C) The eigenvalue for the first PC is removed and the
eigenvalues corresponding to the remaining PCs are close to 0 after the forth
PC. (D) The admixture proportions as estimated with ADMIXTURE.

##  4 Discussion

We have developed a novel approach to assess model fit of PCA and the
admixture model based on structure of the residual correlation matrix. We have
shown that it performs well for simulated and real data, using a suit of
different PCA methods, commonly used in the literature, and the ADMIXTURE
software to estimate model parameters. By assessing the residual correlation
structure visually, one is able to detect model misfit and violation of
modelling assumptions.

The model fit is assessed by comparing visually two matrices of residual
correlation coefficients. The theoretical and practical advantage of our
approach lie in three aspects. First, our approach is computationally simple
and fast. Calculation of the two residual correlation matrices and their
difference is computationally inexpensive. Secondly, our approach provides a
unified approach to model fitting based on PCA and clustering methods (like
ADMIXTURE). In particular, it provides simple means to assess the adequacy of
the chosen number of top principal components to describe the structure of the
data. Assessing the adequacy by plotting the principal components against each
other might lead to false confidence. In contrast, our approach exposes model
misfit by plotting the difference between two matrices of the residual
correlation coefficients. Thirdly, it comes with theoretical guarantees in
some cases. These guarantees are further back up by simulations in cases, we
cannot provide theoretical validity. Finally, our approach might be adapted to
work on NGS data without estimating genotypes first, but working directly on
genotype likelihoods.

## Data availability

The data sets used in this study are all publicly available, including
simulated and real data. Information about the R code used to analyze and
simulate data is available at https://github.com/Ginwaitthreebody/evalPCA. The
variant calls for the 1000 Genomes Project data used are publicly available at
ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/.

## Acknowledgements

The authors are supported by the Independent Research Fund Denmark (grant
number: 8021-00360B) and the University of Copenhagen through the Data+
initiative. SL acknowledges the financial support from the funding agency of
China Scholarship Council. GGE and AA are supported by the Independent
Research Fund Denmark (grant numbers: 8049-00098B and DFF-0135-00211B
respectively).

## Appendix A

We first state the expectation and covariance matrix of $G_{s\star}$ and
$\Pi_{s\star}$, respectively, under the given distributional assumptions,

$\displaystyle\operatorname{\mathbb{E}}[\Pi_{s\star}]$
$\displaystyle=\mu^{T}Q,\qquad\operatorname{cov}(\Pi_{s\star})=Q^{T}\Sigma Q,$
$\displaystyle\operatorname{\mathbb{E}}[G_{s\star}]$
$\displaystyle=2\operatorname{\mathbb{E}}[F_{s\star}Q]=2\mu^{T}Q,$
$\displaystyle\operatorname{cov}(G_{s\star})$
$\displaystyle=\operatorname{\mathbb{E}}[\operatorname{cov}(G_{s\star}\mid\Pi)]+4\operatorname{cov}(\Pi_{s\star})=D+4Q^{T}\Sigma
Q,$

for $s=1,\ldots,m$, where

$\displaystyle D$
$\displaystyle=2\operatorname{\mathbb{E}}[\operatorname{diag}(\Pi_{s1}(1-\Pi_{s1}),\ldots,\Pi_{sn}(1-\Pi_{sn}))],$

and

$\operatorname{\mathbb{E}}[\Pi_{si}(1-\Pi_{si})]=\mu^{T}Q_{\star
i}-(\mu^{T}Q_{\star i})^{2}-(Q^{T}\Sigma Q)_{ii}.$

The unconditional columns $G_{s\star}$, $s=1,\ldots,m$, of $G$ are independent
random vectors by construction.

The above implies that

$\frac{1}{m}\operatorname{\mathbb{E}}[G^{T}G]=D+4Q^{T}(\Sigma+\mu\mu^{T})Q.$
(4)

Auxiliary results are in appendix B.

######  Lemma 7.

The estimator $\widehat{D}$ is an unbiased estimator of $D$, that is,
$\operatorname{\mathbb{E}}[\widehat{D}]=D$. Furthermore, it holds that
$\widehat{D}\to D$ as $m\to\infty$.

###### Proof.

Conditional on $\Pi_{si}$, using binomiality, we have
$\operatorname{\mathbb{E}}[G_{si}(2-G_{si})\mid\Pi_{si}]=2\Pi_{si}(1-\Pi_{si})$,
and the first result follows. For convergence, note that $G_{si}(2-G_{si})$,
$s=1,\ldots,m$, unconditionally, form a sequence of iid random variables with
finite variance, hence the convergence statement follows from the strong Law
of Large Numbers (Jacod and Protter, 2004). ∎

######  Lemma 8.

The estimator $\widehat{H}=\frac{1}{m}G^{T}G-\widehat{D}$ is an unbiased
estimator of $H=4Q^{T}(\Sigma+\mu\mu^{T})Q$, that is,
$\operatorname{\mathbb{E}}[\widehat{H}]=H$. Furthermore, it holds that
$\widehat{H}\to 4Q^{T}(\Sigma+\mu^{T}\mu)Q$ as $m\to\infty$, and

$\operatorname{\mathbb{E}}\left[\|\widehat{H}-4Q^{T}(\Sigma+\mu\mu^{T})Q\|_{F}^{2}\right]\leq\frac{16n^{2}}{m}.$

###### Proof.

Unbiasedness follows from (4) and Lemma 7. Consider the $(i,j)$-th entry of
$\tfrac{1}{m}G^{T}G$, namely, $\tfrac{1}{m}\sum_{s=1}^{m}G_{si}G_{sj}$. The
sequence $G_{si}G_{sj}$, $s=1,\ldots,m$, is iid with finite variance, hence
$\tfrac{1}{m}G^{T}G$ converges to $\operatorname{\mathbb{E}}[G^{T}G]$ as
$m\to\infty$ by the strong Law of Large Numbers (Jacod and Protter, 2004).
Combined with Lemma 7 gives convergence of $\widehat{H}$ to $H$ as
$m\to\infty$.

It remains to prove the inequality. Define

$\displaystyle
A_{s,ij}=\begin{cases}G_{si}G_{sj}-4(Q^{T}(\Sigma+\mu\mu^{T})Q)_{ij}&\text{ if
}i\neq j,\\\ 2G_{si}(G_{si}-1)-4(Q^{T}(\Sigma+\mu\mu^{T})Q)_{ii}&\text{ if
}i=j.\end{cases}$

Then,

$\displaystyle(\widehat{H}_{ij}-\operatorname{\mathbb{E}}[\widehat{H}_{ij}])^{2}$
$\displaystyle=\left(\frac{1}{m}\sum_{s=1}^{m}A_{s,ij}\right)^{2}=\frac{1}{m^{2}}\sum_{s=1}^{m}\sum_{t=1}^{m}A_{s,ij}A_{t,ij},$
$\displaystyle\|\widehat{H}-\operatorname{\mathbb{E}}[\widehat{H}]\|_{F}^{2}$
$\displaystyle=\frac{1}{m^{2}}\sum_{i=1}^{n}\sum_{j=1}^{n}\sum_{s=1}^{m}\sum_{t=1}^{m}A_{s,ij}A_{t,ij}.$

Using $\operatorname{\mathbb{E}}[A_{s,ij}]=0$, independence of $A_{s,ij}$ and
$A_{t,ij}$ for $s\neq t$, and $|A_{s,ij}|\leq 4$, we have

$\displaystyle\operatorname{\mathbb{E}}[\|\widehat{H}-\operatorname{\mathbb{E}}[\widehat{H}]\|_{F}^{2}]$
$\displaystyle=\frac{1}{m^{2}}\sum_{i=1}^{n}\sum_{j=1}^{n}\sum_{s=1}^{m}\operatorname{\mathbb{E}}[A_{s,ij}^{2}]\leq\frac{1}{m^{2}}16mn^{2}=\frac{16n^{2}}{m},$

which proves the claim. ∎

The convergence result is also in Chen and Storey (2015, theorem 2). The
second part provides the rate of convergence of $\widehat{H}$ in the
$L^{2}$-norm. Convergence is contingent on large $m$, rather than large $n$,
and requires $m$ to increase at least like the square of $n$.

Proof of Theorem 1. Since $\widehat{P}_{k^{\prime}}$ is assumed to be an
orthogonal projection, that is,
$\widehat{P}_{k^{\prime}}^{2}=\widehat{P}_{k^{\prime}}$ and
$\widehat{P}_{k^{\prime}}^{T}=\widehat{P}_{k^{\prime}}$, then also the limit
is an orthogonal projection, $P_{k^{\prime}}^{2}=P_{k^{\prime}}$ and
$P_{k^{\prime}}^{T}=P_{k^{\prime}}$.

Consider the empirical covariance $\widehat{B}$. Define the variables
$T_{k^{\prime}}=G(I-P_{k^{\prime}})$ with $\widehat{P}_{k^{\prime}}$ replaced
by $P_{k^{\prime}}$, and the empirical covariance

$\displaystyle\widetilde{B}_{ij}$
$\displaystyle=\frac{1}{m-1}\sum_{s=1}^{m}(T_{k^{\prime},si}T_{k^{\prime},sj}-\hbox{\vbox{\hrule
height=0.5pt\kern 1.29167pt\hbox{\kern-1.00006pt$T$\kern
0.0pt}}}_{k^{\prime},i}\hbox{\vbox{\hrule height=0.5pt\kern
1.29167pt\hbox{\kern-1.00006pt$T$\kern 0.0pt}}}_{k^{\prime},j})$
$\displaystyle=\frac{1}{m-1}\sum_{s=1}^{m}T_{k^{\prime},si}T_{k^{\prime},sj}-\frac{m}{m-1}\hbox{\vbox{\hrule
height=0.5pt\kern 1.29167pt\hbox{\kern-1.00006pt$T$\kern
0.0pt}}}_{k^{\prime},i}\hbox{\vbox{\hrule height=0.5pt\kern
1.29167pt\hbox{\kern-1.00006pt$T$\kern 0.0pt}}}_{k^{\prime},j},$

defined similarly to $\widehat{B}_{ij}$, with $\hbox{\vbox{\hrule
height=0.5pt\kern 1.29167pt\hbox{\kern-1.00006pt$T$\kern
0.0pt}}}_{k^{\prime},i}=\frac{1}{m}\sum_{s=1}^{m}T_{k^{\prime},si}.$ The
sequences $T_{k^{\prime},si}T_{k^{\prime},sj}$, $s=1,2,\ldots$, and
$T_{k^{\prime},si}$, $s=1,2,\ldots$, are iid random variables, by the
distributional assumptions on $G$. Furthermore, since $P_{k^{\prime}}$ is an
orthogonal projection, then $\|I-P_{k^{\prime}}\|^{2}_{F}\leq n$ is bounded
(Lemma 12). Therefore, also $T_{k^{\prime},si}$ is bounded uniformly in $s,i$
by $2\sqrt{n}\leq 2n$.

Using boundedness, independence and the strong Law of Large Numbers (Jacod and
Protter, 2004),

$\displaystyle\widetilde{B}_{ij}$
$\displaystyle\to\operatorname{\mathbb{E}}[T_{k^{\prime},1i}T_{k^{\prime},1j}]-\operatorname{\mathbb{E}}[T_{k^{\prime},1i}]\operatorname{\mathbb{E}}[T_{k^{\prime},1j}]=\operatorname{cov}(T_{k^{\prime},1i},T_{k^{\prime},1j}),$
(5)

for $m\to\infty$, and
$\operatorname{cov}(T_{k^{\prime},1i},T_{k^{\prime},1j})=(I-P_{k^{\prime}})(D+4Q^{T}\Sigma
Q)(I-P_{k^{\prime}}).$ The latter equality follows from (4).

Consider
$R=G(I-\widehat{P}_{k^{\prime}})=G(I-P_{k^{\prime}})+G(P_{k^{\prime}}-\widehat{P}_{k^{\prime}})=T+G(P_{k^{\prime}}-\widehat{P}_{k^{\prime}})$.
Hence,

$\displaystyle{\large|}\ {\hbox{\vbox{\hrule height=0.5pt\kern
1.29167pt\hbox{\kern-1.00006pt$R$\kern
0.0pt}}}}_{k^{\prime},i}-{\hbox{\vbox{\hrule height=0.5pt\kern
1.29167pt\hbox{\kern-1.00006pt$T$\kern 0.0pt}}}}_{k^{\prime},i}{\large|}$
$\displaystyle\leq\frac{1}{m}\sum_{s=1}^{m}\sum_{i^{\prime}=1}^{n}\sum_{j^{\prime}=1}^{n}2|(P_{k^{\prime}}-\widehat{P}_{k^{\prime}})_{i^{\prime}j^{\prime}}|$
$\displaystyle=2\sum_{i^{\prime}=1}^{n}\sum_{j^{\prime}=1}^{n}|(P_{k^{\prime}}-\widehat{P}_{k^{\prime}})_{i^{\prime}j^{\prime}}|\to
0,$

as $m\to\infty$ by assumption of the theorem. It follows that
${\hbox{\vbox{\hrule height=0.5pt\kern 1.29167pt\hbox{\kern-1.00006pt$R$\kern
0.0pt}}}}_{k^{\prime},i}$ converges to
$\operatorname{\mathbb{E}}[T_{k^{\prime},1i}]$ as $m\to\infty$. Furthermore,

$\displaystyle\frac{1}{m-1}\sum_{s=1}^{m}R_{k^{\prime},si}R_{k^{\prime},sj}-\frac{1}{m-1}\sum_{s=1}^{m}T_{k^{\prime},si}T_{k^{\prime},sj}$
$\displaystyle=\frac{1}{m-1}\sum_{s=1}^{m}(T_{si}+(G(P_{k^{\prime}}-\widehat{P}_{k^{\prime}}))_{si})(T_{sj}+(G(P_{k^{\prime}}-\widehat{P}_{k^{\prime}}))_{sj})$
$\displaystyle\quad-\frac{1}{m-1}\sum_{s=1}^{m}T_{k^{\prime},si}T_{k^{\prime},sj}$
$\displaystyle=\frac{1}{m-1}\sum_{s=1}^{m}T_{si}(G(P_{k^{\prime}}-\widehat{P}_{k^{\prime}}))_{sj}+\frac{1}{m-1}\sum_{s=1}^{m}(G(P_{k^{\prime}}-\widehat{P}_{k^{\prime}}))_{si}T_{sj}$
$\displaystyle\quad+\frac{1}{m-1}\sum_{s=1}^{m}(G(P_{k^{\prime}}-\widehat{P}_{k^{\prime}}))_{si}(G(P_{k^{\prime}}-\widehat{P}_{k^{\prime}}))_{sj}.$

The absolute value of the first term in the last line above is bounded by

$\frac{4nm}{m-1}\sum_{j^{\prime}=1}^{n}|(P_{k^{\prime}}-\widehat{P}_{k^{\prime}})_{j^{\prime}j}|,$

and similarly for the second term. The third is bounded by

$\frac{4m}{m-1}\sum_{i^{\prime}=1}^{n}\sum_{j^{\prime}=1}^{n}|(P_{k^{\prime}}-\widehat{P}_{k^{\prime}})_{i^{\prime}i}||(P_{k^{\prime}}-\widehat{P}_{k^{\prime}})_{j^{\prime}j}|.$

All three terms converge to zero as $m\to\infty$, hence we conclude from (5)
that
$\widehat{B}_{ij}\to\operatorname{cov}(T_{k^{\prime},1i},T_{k^{\prime},1j})$
as $m\to\infty$.

The result for the estimated covariance $\widehat{C}$ follows from convergence
of $\widehat{D}$ and by assumption of the theorem. The remaining part follows
from the convergence of $\widehat{B}$ and $\widehat{C}$. Note that $QP=Q$,
hence the second equation holds. The last statement of the theorem follows
directly.

Proof of Theorem 2. Consider $T_{k}=G(I-P_{k})=G(1-P)$, where
$P=Q^{T}(QQ^{T})^{-1}Q$ is the projection onto the row space of $Q$. Then,
$T_{k}$ contains the residuals under multiple regression of the $m$ rows of
$G$ on the $k$ rows of $Q$ (Box et al., 2005). Since $e$ is in the row space
of $Q$, then the sum of the residuals is zero for each $s=1,\ldots,m$:
$\sum_{i=1}^{n}T_{k,si}=0$ (the assumption that $e$ is in the row space is
equivalent to having an intercept in the regression model) (Box et al., 2005).
We have, for $s=1,\ldots,m$,

$\displaystyle 0$
$\displaystyle=\operatorname{var}\left(\sum_{i=1}^{n}T_{k,si}\right)=\sum_{i=1}^{n}\operatorname{var}(T_{k,si})+\sum_{i=1}^{n}\sum_{j=1,i\not=j}^{n}\operatorname{cov}(T_{k,si},T_{k,sj})$
$\displaystyle=\sum_{i=1}^{n}\operatorname{var}(T_{k,1i})+\sum_{i=1}^{n}\sum_{j=1,i\not=j}^{n}\operatorname{cov}(T_{k,1i},T_{k,1j}),$

since the distribution of $T_{k,si}$ is independent of $s$. From the proof of
Theorem 4, it follows that $\widehat{B}$ converges to
$\operatorname{cov}(T_{k,1\star})$ as $m\to\infty$. Hence,

$\displaystyle\sum_{i=1}^{n}\widehat{B}_{ii}+\sum_{i=1}^{n}\sum_{j=1,i\not=j}^{n}\widehat{B}_{ij}$
$\displaystyle\to\ 0,\quad\text{as}\quad m\to\infty,$

and the desired result follows by rearrangement.

If $Q$ takes the given form, then the residuals under multiple regression are
independent between compartments, as the projection is

$P=\begin{pmatrix}P_{1}&0&\cdots&0\\\ 0&P_{2}&\cdots&0\\\
\vdots&\vdots&\ddots&\vdots\\\ 0&0&\cdots&P_{r}\end{pmatrix},$

where $P_{\ell}=Q_{\ell}^{T}(Q_{\ell}Q_{\ell}^{T})^{-1}Q_{\ell}$ has dimension
$n_{\ell}\times n_{\ell}$. It follows that the computation above holds for
each compartment. Finally, if $Q_{\ell}=(1\ldots 1)$, then the distribution of
the random vector $T_{k,1\star}$ is exchangeable, resulting in

$\displaystyle 0$
$\displaystyle=\operatorname{var}\left(\sum_{i=1}^{n_{\ell}}T_{k,1i}\right)=\sum_{i=1}^{n_{\ell}}\operatorname{var}(T_{k,1i})+\sum_{i=1}^{n_{\ell}}\sum_{j=1,i\not=j}^{n_{\ell}}\operatorname{cov}(T_{k,1i},T_{k,1j})$
$\displaystyle=n_{\ell}\operatorname{var}(T_{k,11})+n_{\ell}(n_{\ell}-1)\operatorname{cov}(T_{k,11},T_{k,12})$

assuming the individuals in the $\ell$-th compartment are numbered $1$ to
$n_{\ell}$. Rearranging terms and substituting $\widehat{b}_{ij}$ for the
moments of $T_{k,i\star}$ yields the desired result.

Proof of Theorem 3. Consider $T_{k}=G(I-P_{k})=G(1-P)$, where
$P=Q^{T}(QQ^{T})^{-1}Q$ is the projection onto the row space of $Q$. If
$Q_{1}=(1\ldots 1)$, then the distribution of the random variables
$T_{k,11},\ldots,T_{k,1n_{1}}$ are exchangeable, resulting in

$\displaystyle 0$
$\displaystyle\leq\operatorname{var}\left(\sum_{i=1}^{n_{1}}T_{k,1i}\right)=\sum_{i=1}^{n_{1}}\operatorname{var}(T_{k,1i})+\sum_{i=1}^{n_{1}}\sum_{j=1,i\not=j}^{n_{1}}\operatorname{cov}(T_{k,1i},T_{k,1j})$
$\displaystyle=n_{1}\operatorname{var}(T_{k,11})+n_{1}(n_{1}-1)\operatorname{cov}(T_{k,11},T_{k,12}).$

Rearranging terms and substituting $\widehat{b}_{ij}$ for the moments of
$T_{k,i\star}$ yields the desired result.

Proof of Theorem 4. The convergence statement of the theorem is a special case
of Theorem 9 in Appendix B. Take $A_{m}=\widehat{H}$ (that depends on the
number of SNPs $m$, and the particular realization), $A=H$, and $k=k^{\prime}$
in the theorem ($k$ is used as a generic index in Theorem 9). Then,
$E_{k}E_{k}^{T}=P_{k^{\prime}}$ and
$F_{m,k^{\prime}}F_{m,k}^{T}=\widehat{P}_{k^{\prime}}$, and the conclusion of
Theorem 4 holds. Convergence in Frobenius norm is equivalent to pointwise
convergence (as $n$ is fixed) $\widehat{P}_{k^{\prime}}\to P_{k^{\prime}}$ as
$m\to\infty$ by definition.

If $\Sigma+\mu\mu^{T}$ is positive definite, then it has rank $k$. As
$\text{rank}(Q)=k$ by assumption, it follows from Lemma 13 that
$\text{rank}(H)=k$. Consequently, there are $k$ positive eigenvalues of $H$
and $\lambda_{k+1}=0$, and the eigenvalue condition holds. Conversely, assume
the eigenvalue condition holds. By definition $\text{rank}(H)\leq k$. As
$\lambda_{k}>\lambda_{k+1}\geq 0$ by assumption, then also $\text{rank}(H)\geq
k$ and we conclude $\text{rank}(H)=k$. It follows that the rank of
$\Sigma+\mu\mu^{T}$ is $k$; consequently, it is positive definite.

If $k^{\prime}=k=n$, then $P_{k}=V_{k}V_{k}^{T}=I$ and $P=I$ (as $k=n$), and
$P_{k}=P$. So assume $k^{\prime}=k<n$. Since the eigenvalue condition is
fulfilled, then from the above, we have
$\text{rank}(Q^{T}(\Sigma+\mu\mu^{T}))=k$, and Lemma 13 yields that the row
space of $H$ and $Q$ agree. Similarly, we have
$H=V_{k}\operatorname{diag}(\lambda_{1},\ldots,\lambda_{k})V^{T}_{k}$ and
Lemma 13 yields that the row space of $H$ and $V_{k}^{T}$ agree. This implies
the row space of $Q$ and $V_{k}^{T}$ agree. Consequently,
$P_{k}=Q^{T}(QQ^{T})^{-1}Q=P$, and the statement holds.

Proof of Theorem 5. It follows trivially that $e$ is an eigenvector of $H_{1}$
with eigenvalue $0$. If $D$ has all entries positive, then it is positive
definite and $D+4Q^{T}(\Sigma+\mu\mu^{T})Q$ is also positive definite, hence
has rank $n$. It follows from Lemma 13 that $H_{1}$ has rank $n-1$, hence
$\lambda_{n-1}>0$.

Similarly to the proof of Lemma 8 in Appendix B, one can show
$\operatorname{\mathbb{E}}[\widehat{H}_{1}]=H_{1}$ and $\widehat{H}_{1}\to
H_{1}$ as $m\to\infty$, where $H_{1}$ denotes the right hand side of (3). The
remaining part of the theorem is proven similarly to Theorem 4.

Proof of Theorem 6. Note that $e$ is an eigenvector of $H_{1}$ with eigenvalue
$0$. Consider an eigenvector $v$ of $H_{1}$, orthogonal to $e$ with eigenvalue
$\lambda$. Then, the following two equations are equivalent,

$\displaystyle\left(I-\frac{1}{n}E\right)(D+4Q^{T}(\Sigma+\mu\mu^{T})Q)\left(I-\frac{1}{n}E\right)v$
$\displaystyle=\lambda v,$ $\displaystyle
4\left(I-\frac{1}{n}E\right)Q^{T}(\Sigma+\mu\mu^{T})Q\left(I-\frac{1}{n}E\right)v$
$\displaystyle=(\lambda-d)v,\ $ (6)

where it is used that $D=dI$ and $v\perp e$. It shows that $v$ is an
eigenvector of
$K=4(I-\tfrac{1}{n}E)Q^{T}(\Sigma+\mu\mu^{T})Q(I-\tfrac{1}{n}E)$ with
eigenvalue $\mu=\lambda-d$. Since $Q$ has rank $k$ and the vector $e$ is in
the space spanned by the rows of $Q$, then $Q(I-\tfrac{1}{n}E)$ has rank
$k-1$. It follows that there are at most $k-1$ positive eigenvalues of $K$,
that is, at most $k-1$ eigenvalues of $H_{1}$ such that $\lambda>d$.
Furthermore, there are precisely $k-1$ positive eigenvalues, provided
$\Sigma+\mu\mu^{T}$ is positive definite (Lemma 13). The remaining eigenvalues
of $K$ are zero, that is, the corresponding eigenvalues of $H_{1}$ are
$\lambda=d$.

Assume $\Sigma+\mu\mu^{T}$ is positive definite, then by the above argument
there precisely are $k-1$ eigenvalues of $H_{1}$ such that $\lambda>d$ with
corresponding orthogonal eigenvectors $v_{1},\ldots,v_{k-1}$. It follows from
(6) that $v_{1},\ldots,v_{k-1}$ are in the space spanned by the rows of
$Q(I-\tfrac{1}{n}E)$, hence the eigenvectors are in the space spanned by the
rows of $Q$. By assumption $e$ is also in that row span. Hence,
$v_{1},\ldots,v_{k-1},e$ forms an orthogonal basis of the row span of $Q$, as
$Q$ has rank $k$. Thus, $P_{k}=P$.

## Appendix B

######  Theorem 9.

Let $A_{m}$ be a sequence of symmetric $n\operatorname{\times}n$-matrices that
converges to a symmetric $n\operatorname{\times}n$-matrix $A$ in the Frobenius
norm, that is $\|A_{m}-A\|_{F}\to 0$, as $m\to\infty$. Let
$\lambda_{1}\geq\ldots\geq\lambda_{n}$ be the eigenvalues of $A$ (with
multiplicity, and not necessarily non-negative). Let $k\leq n$ be given and
assume either $k=n$ or $\lambda_{k}>\lambda_{k+1}$. Furthermore, let
$e_{1},\ldots,e_{k}$ be orthogonal eigenvectors corresponding to the
eigenvalues $\lambda_{1},\ldots,\lambda_{k}$, respectively, and let
$f_{m,1},\ldots,f_{m,k}$ be orthogonal eigenvectors corresponding to the $k$
largest eigenvalues of $A_{m}$ (with multiplicity). Then, the orthogonal
projection onto the span of $f_{m,1},\ldots,f_{m,k}$ converges to the
orthogonal projection onto the span of $e_{1},\ldots,e_{k}$ in the Frobenius
norm. That is, define $E_{k}=(e_{1},\ldots,e_{k})$ and
$F_{m,k}=(f_{m,1},\ldots,f_{m,k})$, then
$\|F_{m,k}F_{m,k}^{T}-E_{k}E_{k}^{T}\|_{F}\to 0$ as $m\to\infty$.

###### Proof.

If $k=n$, then $E_{n}E_{n}^{T}=I$ and $F_{m,n}F_{m,n}^{T}=I$, and the
statement is trivial. Hence, assume $k<n$. Let $e_{1},\ldots,e_{n}$ be
eigenvectors of $A$ corresponding to eigenvalues
$\lambda_{1},\ldots,\lambda_{n}$, respectively. Let $f_{m,1},\ldots,f_{m,n}$
be the eigenvectors of $A_{m}$ corresponding to the eigenvalues
$\mu_{m,1}\geq\ldots\geq\mu_{m,n}$. All eigenvectors can be asssumed to be
orthonormal.

As $\|A-A_{m}\|_{F}^{2}\to 0$ for $m\to\infty$, then every entry of $A_{m}$
converges to the corresponding entry of $A$. Consequently, the characteristic
function of $A_{m}$ converges to that of $A$, and the eigenvalues of $A_{m}$
converges to those of $A$, that is, $\mu_{m,j}\to\lambda_{j}$ for
$j=1,\ldots,n$, and $m\to\infty$. Let $T_{m}$ be such that
$E_{n}=F_{m,n}T_{m}$. As $E_{n}$ and $F_{m,n}$ are orthogonal matrices, hence
also $T_{m}$ is orthogonal. Applying Lemma 10 in the first and third line
gives

$\displaystyle\|A-A_{m}\|_{F}^{2}$ $\displaystyle=\|AE-
A_{m}E_{n}\|_{F}^{2}=\|E\operatorname{diag}(\lambda_{1},\ldots,\lambda_{n})-A_{m}F_{m,n}T_{m}\|_{F}^{2}$
$\displaystyle=\|F_{m,n}T_{m}\operatorname{diag}(\lambda_{1},\ldots,\lambda_{n})-F_{m,n}\operatorname{diag}(\mu_{m,1},\ldots,\mu_{m,n})T_{m}\|_{F}^{2}$
$\displaystyle=\|T_{m}\operatorname{diag}(\lambda_{1},\ldots,\lambda_{n})-\operatorname{diag}(\mu_{m,1},\ldots,\mu_{m,n})T_{m}\|_{F}^{2}$
$\displaystyle=\sum_{i=1}^{n}\sum_{j=1}^{n}(\lambda_{j}T_{m,ij}-\mu_{m,i}T_{m,ij})^{2}=\sum_{i=1}^{n}\sum_{j=1}^{n}T_{m,ij}^{2}(\lambda_{j}-\mu_{m,i})^{2}.$

By assumption, $\lambda_{k}>\lambda_{k+1}$. Hence, by convergence of
eigenvalues, for $j\leq k$, $i\geq k+1$, or $j\geq k+1$, $i\leq k$, we have
$T_{m,ij}\to 0$ as $m\to\infty$.

Furthermore,

$\displaystyle E_{k}E_{k}^{T}-F_{m,k}F_{m,k}^{T}$
$\displaystyle=\sum_{\ell=1}^{k}\left(e_{\ell}e_{\ell}^{T}-f_{m,\ell}f_{m,\ell}^{T}\right)$
$\displaystyle=\sum_{\ell=1}^{k}\Big{(}\Big{(}\sum_{a=1}^{n}f_{m,a}T_{m,a\ell}\Big{)}\Big{(}\sum_{a=1}^{n}f_{m,a}T_{m,a\ell}\Big{)}^{\\!\\!T}-f_{m,\ell}f_{m,\ell}^{T}\Big{)}$
$\displaystyle=\sum_{\ell=1}^{k}\left(\sum_{a=1}^{n}\sum_{b=1}^{n}T_{m,a\ell}T_{m,b\ell}f_{m,a}f_{m,b}^{T}-f_{m,\ell}f_{m,\ell}^{T}\right)$
$\displaystyle=\sum_{a=1}^{n}\sum_{b=1}^{n}\sum_{\ell=1}^{k}T_{m,a\ell}T_{m,b\ell}f_{m,a}f_{m,b}^{T}-\sum_{\ell=1}^{k}f_{m,\ell}f_{m,\ell}^{T}.$
$\displaystyle=\sum_{(a,b)\in\\{1,\ldots,n\\}^{2}\setminus
A_{1,k}}\Big{(}\sum_{i=1}^{k}T_{m,ai}T_{m,bi}\Big{)}f_{m,a}f_{m,b}^{T}$
$\displaystyle\quad+\sum_{(a,a)\in
A_{1,k}}\Big{(}\sum_{i=1}^{k}T_{m,ai}T_{m,ai}-1\Big{)}f_{m,a}f_{m,a}^{T},$

where $A_{i,j}=\\{(a,a):i\leq a\leq j\\}$.

From Lemma 11, we have $f_{m,a}f_{m,b}^{T}\perp f_{m,c}f_{m,d}^{T}$ for
$(a,b)\neq(c,d)$ in the Frobenius inner product. Moreover,
$\|f_{m,a}f_{m,b}^{T}\|_{F}=1$ for all $a,b$. Hence,

$\displaystyle\|E_{k}E_{k}^{T}-F_{m,k}F_{m,k}^{T}\|_{F}^{2}$
$\displaystyle=\sum_{(a,b)\in\\{1,\ldots,n\\}^{2}\setminus
A_{1,k}}\Big{(}\sum_{i=1}^{k}T_{m,ai}T_{m,bi}\Big{)}^{\\!\\!2}+\sum_{(a,a)\in
A_{1,k}}\Big{(}\sum_{i=1}^{k}T_{m,ai}T_{m,ai}-1\Big{)}^{\\!\\!2}.$
$\displaystyle=\sum_{(a,b)\in\\{1,\ldots,n\\}^{2}\setminus
A_{1,n}}\Big{(}\sum_{i=1}^{k}T_{m,ai}T_{m,bi}\Big{)}^{\\!\\!2}+\sum_{(a,a)\in
A_{k+1,n}}\Big{(}\sum_{i=1}^{k}T_{m,ai}T_{m,ai}\Big{)}^{\\!\\!2}$
$\displaystyle\quad+\sum_{(a,a)\in
A_{1,k}}\Big{(}\sum_{i=1}^{k}T_{m,ai}T_{m,ai}-1\Big{)}^{\\!\\!2}.$ (7)

As noted above, $T_{m,ij}\to 0$ as $m\to\infty$ for $j\leq k$, $i\geq k+1$, or
$j\geq k+1$, $i\leq k$. Using this and orthogonality of $T_{m}$ gives

$\sum_{i=1}^{k}T_{m,ai}T_{m,bi}=\sum_{i=1}^{n}T_{m,ai}T_{m,bi}-\sum_{i=k+1}^{n}T_{m,ai}T_{m,bi}\to\begin{cases}0&\text{
if }a\neq b,\\\ 1&\text{ if }a=b,\end{cases}.$

Inserting into (7) results in
$\|E_{k}E_{k}^{T}-F_{m,k}F_{m,k}^{T}\|_{F}^{2}\to 0,$ as $m\to\infty$. ∎

######  Lemma 10.

Let $A$ be an $a\operatorname{\times}b$ matrix. Let $U$ be a
$b\operatorname{\times}b$ orthogonal matrix and $V$ an
$a\operatorname{\times}a$ orthogonal matrix. Then,

$\displaystyle\|A\|_{F}=\|VA\|_{F}=\|AU\|_{F}=\|VAU\|_{F}.$

###### Proof.

See Golub and Loan (2013). ∎

######  Lemma 11.

Let $w,x,y,z\in\operatorname{\mathbb{R}}^{b}$. For
$a\operatorname{\times}b$-matrices $A$ and $B$, let $\left\langle
A,B\right\rangle_{F}=\sum_{i=1}^{a}\sum_{j=1}^{b}A_{ij}B_{ij}$ be the
Frobenius inner product of $A$ and $B$, and let
$\left\langle\operatorname{\,\cdot\,},\operatorname{\,\cdot\,}\right\rangle$
be the standard inner product on $\operatorname{\mathbb{R}}^{b}$. Then,
$\left\langle wx^{T},yz^{T}\right\rangle_{F}=\left\langle
w,y\right\rangle\left\langle x,z\right\rangle$. In particular,
$\|wx^{T}\|_{F}=\|w\|_{2}\|x\|_{2}$ and $wx^{T}\perp yz^{T}$ if $w\perp y$ or
$x\perp z$.

###### Proof.

Note that

$\displaystyle\left\langle
wx^{T},yz^{T}\right\rangle=\sum_{i=1}^{b}\sum_{j=1}^{b}w_{i}x_{j}y_{i}z_{j}=\sum_{i=1}^{b}w_{i}y_{i}\sum_{j=1}^{b}x_{j}z_{j}=\left\langle
w,y\right\rangle\left\langle x,z\right\rangle.$

Hence, if either $w\perp y$ or $x\perp z$, then $wx^{T}\perp yz^{T}$, and
$\|wx^{T}\|_{F}^{2}=\left\langle wx^{T},wx^{T}\right\rangle=\left\langle
w,w\right\rangle\left\langle x,x\right\rangle=\|w\|_{2}^{2}\|x\|_{2}^{2}$,
such that $\|wx^{T}\|_{F}=\|w\|_{2}\|x\|_{2}$. ∎

######  Lemma 12.

Let $v_{1},\ldots,v_{\ell}$ be linearly independent vectors. An orthogonal
projection matrix on $span(v_{1},\ldots,v_{\ell})$ has Frobenius norm
$\sqrt{\ell}$.

###### Proof.

We may assume that $v_{1},\ldots,v_{\ell}$ are orthonormal. Then, we can write
the projection matrix as $P=v_{1}v_{1}^{T}+\ldots+v_{\ell}v_{\ell}^{T}$. By
Lemma 11, $v_{i}v_{i}^{T}\perp_{F}v_{j}v_{j}^{T}$ for $i\neq j$. So, again by
Lemma 11,
$\|P\|_{F}^{2}=\|v_{1}v_{1}^{T}\|_{F}^{2}+\ldots+\|v_{\ell}v_{\ell}^{T}\|_{F}^{2}=\ell$.
∎

######  Lemma 13.

Let $A$ be an $a\times b$ matrix and $B$ an $b\times c$ matrix, both of rank
$b$, such that $a,c\geq b$. Let $C=AB.$ Then, $C$ is of rank $b$, and the row
space of $C$ coincides with the row space of $B$.

###### Proof.

First we show that $\text{rank}(C)=b$. Note that $A$ has $b$ linearly
independent rows $1\leq i_{1}<\ldots<i_{b}\leq b$, and $B$ has $b$ linearly
independent columns $1\leq j_{1}<\ldots<j_{b}\leq b$. Let $\widetilde{A}$ and
$\widetilde{B}$ be the $b\times b$ matrices with
$\widetilde{A}_{cd}=A_{i_{c}d}$ and $\widetilde{B}_{cd}=B_{cj_{d}}$. Then
$\widetilde{A}$ and $\widetilde{B}$ are invertible matrices. Hence, also
$\widetilde{C}=\widetilde{A}\widetilde{B}=(C_{i_{a}j_{b}})_{a,b}$ is
invertible and has rank $k$. It follows that $C$ has rank $k$. As
$\widetilde{A}$ is invertible, then the span of the rows of $\widetilde{A}B$
is equal to the span of the rows of $B$. That is, the span of the rows of $AB$
is equal to the span of the rows of $B$. ∎

## References

  * Alexander and Lange (2011) Alexander DH, Lange K. 2011. Enhancement of the admixture algorithm for individual ancestry estimation. BMC Bioinformatics. 12:246.
  * Alexander et al. (2009) Alexander DH, Novembre J, Lange K. 2009. Fast model-based estimation of ancestry in unrelated individuals. Genome Res. 19:1655–1664.
  * Auton et al. (2015) Auton A, Brooks LD, Durbin RM, Garrison EP, Kang HM, Korbel JO, Marchini JL, McCarthy S, McVean GA, Abecasis GR et al. 2015. A global reference for human genetic variation. Nature. 526:68–74.
  * Balding and Nichols (1995) Balding DJ, Nichols RA. 1995. A method for quantifying differentiation between populations at multi-allelic loci and its implications for investigating identity and paternity. Genetica. 96:3–12.
  * Box et al. (2005) Box G, Hunter J, Hunter W. 2005. Statistics for Experimenters: Design, Innovation, and Discovery. Wiley Series in Probability and Statistics. Wiley.
  * Cabreros and Storey (2019) Cabreros I, Storey J. 2019. A Likelihood-Free Estimator of Population Structure Bridging Admixture Models and Principal Components Analysis. Genetics. 212:1009–1029.
  * Chen and Storey (2015) Chen X, Storey J. 2015. Consistent estimation of low-dimensional latent structure in high-dimensional data.
  * Conomos et al. (2016) Conomos M, Reiner A, Weir B, Thornton T. 2016. Model-free estimation of recent genetic relatedness. Am J Hum Genet. 98:127–148.
  * Engelhardt and Stephens (2010) Engelhardt B, Stephens M. 2010. Analysis of population structure: a unifying framework and novel methods based on sparse factor analysis. PLoS Genetics. 6\.
  * Evanno et al. (2005) Evanno G, Regaut S, Goudet J. 2005. Detecting the number of clusters of individuals using the software structure: A simulation study. Mol Ecol. 14:2622–2620.
  * Garcia-Erill and Albrechtsen (2020) Garcia-Erill G, Albrechtsen A. 2020. Evaluation of model fit of inferred admixture proportions. Molecular Ecology Resources. 20:936–949.
  * Golub and Loan (2013) Golub GH, Loan CF. 2013. Matrix Computations. Johns Hopkins Studies in Mathematical Sciences. JHU Press.
  * Jacod and Protter (2004) Jacod J, Protter P. 2004. Probability Essentials. Universitext. Springer.
  * Janes et al. (2017) Janes JK, Miller JM, Dupuis JR, Malenfant RM, Gorrell JC, Cullingham CI, Andrew RL. 2017. The $k02$ conundrum. Mol Ecol. 26:3594–3602.
  * Jolliffe (2002) Jolliffe IT. 2002. Principle Component Analysis (2nd Ed.). Springer Series in Statistics. Springer.
  * Jolliffe and Cadima (2016) Jolliffe T, Cadima J. 2016. Principal component analysis: a review and recent developments. Phil. Trans. R. Soc. A. 374:0150202.
  * Lawson et al. (2018a) Lawson D, van Dorp L, Falush D. 2018a. A tutorial on how not to over-interpret structure and admixture bar plots. Nature Communications. 9.
  * Lawson et al. (2018b) Lawson DJ, van Dorp L, Falush D. 2018b. A tutorial on how not to over-interpret structure and admisture bar plots. Nat Comm. 19:3258.
  * Lazaridis et al. (2014) Lazaridis I, Patterson N, Mittnik A, Renaud G, Mallick S, Kirsanow K, Sudmant PH, Schraiber JG, Castellano S, Lipson M et al. 2014. Ancient human genomes suggest three ancestral populations for present-day Europeans. Nature. 513:409–413.
  * Meisner et al. (2021) Meisner J, Liu S, Huang M, Albrechtsen A. 2021. Large-scale inference of population structure in presence of missingness using PCA. Bioinformatics. 37:1868–1875.
  * Ochoa and Storey (2019) Ochoa A, Storey JD. 2019. $f_{ST}$ and kinship for arbitrary population structures i: Generalized definitions. bioRxiv. .
  * Patterson et al. (2006) Patterson N, Price AL, Reich D. 2006. Population structure and eigenanalysis. PLoS Genetics. 2:e190.
  * Pickrell and Pritchard (2012) Pickrell J, Pritchard J. 2012. Inference of population splits and mixtures from genome-wide allele frequency data. PLOS Genetics. 8:1–17.
  * Pritchard et al. (2000) Pritchard J, Stephens M, Donnelly P. 2000. Inference of population structure using multilocus genotype data. Genetics. 155:945–959.
  * Raj et al. (2014) Raj A, Stephens M, Pritchard J. 2014. Faststructure: Variational inference of populations structure in large snp data sets. Genetics. 197:573–589.
  * Wang (2003) Wang J. 2003. Maximum-likelihood estimation of admixture proportions from genetic data. Genetics. 154:747 –765.
  * Wang (2019) Wang J. 2019. A parsimony estimator of the number of populations froma structure-like analysis. Mol Ecol Res. 19:970 –981.