1. Sparse factor analysis
Let n be the number of individuals in a sample and p be the number of genotypes.
Represent each allele at a locus as a number (e.g., for SNPs from a diploid organism, as in our results above, represent AA as 0, AB as 1 and BB as 2). Our factor
analysis model with K factors can be written as:(1)

Gi,j = µi,j +KX

Λi,k Fk,j + i,j ,k=1

or, equivalently,(2)−1Gi,j ∼ N (µi,j + (ΛF )i,j , ψi,j)

where G is an n×p data matrix, the mean term µi,j is the sum of row- and columnspecific means: µi,j = νi + ξj , Λ is the n × K matrix of factor loadings, F is theK × p matrix of factors, and  is an n × p matrix with each element independently−1
distributed i,j ∼ N (0, ψi,j) and is the product of a row- and a column-specificvariance term ψi,j = θi ηj . We put a gamma prior on the inverse residual variancethat acts as a regularizer: θi ∼ Ga(α, β) and ηj ∼ Ga(κ, τ ), which has mean
αβ (κτ , respectively) and variance αβ 2 (κτ 2 , respectively). In practice, we set
α = κ = 1 and β = 20, τ = 20. This model, with only mean term ξj , is referredpnto as SFAm in the main text; the SFA model is obtained by fixing the vector µat zero. For both SFA and SFAm, as described in the paper, we fix ηj = 1 for
j = 1...p. The ECME algorithm for fitting the general SFA model is describedbelow; the ECME algorithm for fitting SFA is obtained by simply setting µi,j = 0and ηj = 1 throughout. Note that here we have chosen to have column-specific(i.e., SNP-specific) means and row-specific (i.e., individual-specific) variances Ψ.
The other possible options are implemented in the software. In some contexts,
including the population structure problem considered here, it might make sense
to allow more general assumptions, such as variance terms on both the rows andcolumns of the matrix, but we leave their evaluation to future work.
To induce sparsity in the factor loadings Λ, we use an automatic relevance
2determination (ARD) prior (?). Specifically, we assume Λi,k ∼ N (0, σi,k), where2
the matrix Σ = (σi,k )i=1,...,n,k=1,...,K is a parameter that we estimate, together with2
the other parameters, using maximum likelihood. If the estimate of σi,k= 0, thisimplies that Λi,k = 0, thus inducing sparsity.
1

2

Integrating out Λ, the rows of G are conditionally independent given the otherparameters, with:
Gi,· ∼ N (µ, F t Σi F + Ψ−1i ),

(3)

2
2
on the diagonal),) (a diagonal matrix with the K-vector σi,·where Σi = diag(σi,·−1−1 −1−1 −1and Ψi = diag(θi η1 , . . . , θi ηp ). Thus the log marginal likelihood for theparameters µ, F, Σ, Ψ is:L(µ, F, Σ,(4)Ψ; G) := log p(G|µ, F, Σ, Ψ)niX
1h−1 −1tt= −(5)p log(2π) + log |F t Σi F + Ψ−1|+G̃
(FΣF+Ψ)G̃
ii,· ,i,·ii
2i=1where G̃i,j := Gi,j − µi,j .
2. Sparse factor analysis ECME algorithm
We fit this model using an expectation conditional maximization either (ECME)
algorithm (?) to maximize L(µ, F, Σ, Ψ; G). This algorithm is similar to an EM
algorithm, but each maximization step maximizes either the expected log likelihood, or the marginal log likelihood, for a subset of the parameters conditionalon the others. Specifically, the updates to µ, F , and Ψ involve maximizing theexpected log likelihood (with the expectation taken over Λ), whereas the updates
to Σ directly maximize the log marginal likelihood.
To compute the expected log likelihood requires the first and second momentsof the factor loadings Λi,· . The data Gi,· and the loadings Λi,· are jointly normal(as in, e.g., Ghahramani and Hinton (?)):
   t µGi,·F Σi F + Ψ−1F t Σiiµ, F, Σi , Ψi ∼ N,,(6)
0K
Λi,·
Σi F
Σi
where 0K is a K-vector of zeros. Standard results for joint Gaussian distributions
give the conditional expectation for Λi,· :(7)Λi := E [Λi,· |Gi,· , µ, F, Σi , Ψi ] = Ωi G̃i,· ,−1where Ωi = Σi F (F t Σi F + Ψ−1i ) . Similarly, the conditional second moment isgiven by:(8)

Λ2i := E[Λi,· Λti,· |Gi,· , µ, F, Σi , Ψi ] = Σi − Ωi F t Σi + Ωi G̃i,· G̃ti,· Ωti .

The updates for µ, F , and Ψ involve maximizing the expected complete data log
likelihood, Q(µ, F, Σ, Ψ; G) := E[log(p(G|Λ, µ, F, Ψ))|Σ], which from Equation 2is given by:(9)

Q(µ, F, Σ, Ψ; G) = const +nX
i=1Qi (µ, F, Σi , Ψi ; Gi,· )3whereppX+ p(α − 1) log(θi ) ++ p(κ − 1)log(ηj )22j=1
 p θiηj 1X
tt22−.
− ψi
G̃i,j − 2G̃i,j F·,j Λi + F·,j Λi F·,j −2 j=1
τβQi (µ, F, Σi , Ψi ; Gi,· ) =(10)pTaking the derivative of Q(µ, F, Σ, Ψ; G) with respect to νi and setting to 0, we
get the update for νi :(11)p∂Q(F, Σ, Ψ, µ; Gi,· )1Xt
ψi,j −2(Gi,j − µi,j ) + 2Fj,·=Λi = 0∂νi
2 j=1

Pp(12)

ν̂i =tj=1 ηj (Gi,j − ξj − Fj,· Λi )Pp.
j=1 ηjTaking the derivative of Q(µ, F, Σ, Ψ; G) with respect to ξ and setting to 0, we
get the update for ξ:(13)n∂Q(F, Σ, Ψ, µ; G)1X=
Ψi −2(Gi,· − µ) + 2F t Λi = 0∂ξ2 i=1Pn(14)µ̂ =− ξ + F t Λi )i=1 θi (GPi,·n i.
i=1 θiIn these expressions, and in what follows, we are assuming element-wise multiplication when a scalar multiplies a vector or a matrix.
Taking the derivative of Q(µ, F, Σ, Ψ; G) with respect to F·,j and setting to zero,we get the update for F·,j :n∂Q(F, Σ, Ψ, µ; G) X=Ψi (G̃i,j Λi − Λ2i F·,j ) = 0∂F·,ji=1(15)F̂·,j =nX
i=1!−1
θi Λ2inX

θi G̃i,j Λi .

i=1

Taking the derivative of Q(F, Σi , Ψi , µ; Gi,· ) with respect to θi and setting tozero, we get the update for θi :"#p 2 −1X1tt(16) α̂i =ηj G̃ti,j G̃i,j − 2G̃i,j F·,jΛi − F·,jΛ2i F·,j +.
p + 2p(α − 1) j=1
β4

We find the updates for ηj similarly:#"n 2 −1X1tt
θi G̃ti,j G̃i,j − 2G̃i,j F·,jΛi − F·,jΛ2i F·,j +.
(17) η̂j =n + 2n(κ − 1) i=1
τ2To update σi,kwe can use the result from Tipping and Faul (?) to obtain thevalues of Σ that maximize the log marginal likelihood L(µ, F, Σ, Ψ; G) with fixedvalues of µ, F , and Ψ:(18)22σ̂i,k= [(qi,k− si,k )/s2i,k ]+−1−12= Fkt β¬k,iG̃i,· and si,k = Fkt β¬k,iFk , where β¬k,i = (F t Σi,¬k F ) + Ψ−1andwhere qi,ki
2222Σi,¬k = diag(σi,1, ..., σi,k−1, 0, σi,k+1, ..., σi,K). Note that [a]+ = a when a > 0 and= 0 otherwise. This works because, given F , the SFA model (Equation 1) isessentially the sparse regression model considered in Tipping and Faul (?) withF playing the role of the covariates. We implement the simplified version of thisequation from ?. In practice, we recomputed the C matrix (here, the βi matrix) for
each k = 1...K for only the first few iterations of the method, then each iteration
is much faster if we computed βi once for all k = 1...K, with very minimal impact
on the results or convergence time (in the code, the cutoff is set to 5 iterations).
Note that F and Σ are non-identifiable in that multiplying the k th row of F
by a constant c and dividing the k th column of Σ by c2 will not change the likelihood
To deal with this we impose an identifiability
constraint,
Pp (Equation 4).
Pp
112(F−F̄)=1fork=1,...,
K,
whereF̄
=F.
Specifically,k,jk,·k,·k,jj=1
j=1pp
after each iteration we divide every element of Fk,· by its standard deviation ck ,and multiply the k th column of Σ by c2k .
Similarly, the row and column means νj and ξj are non-identifiable in that adding
a constant c to one set of means is equivalent to subtracting c from the other setof means in terms of the matrix product ν1p + ξ1n . Intuitively, one might want tocenter the means for the individuals or samples at zero, allowing genes or SNPs tohave non-centered means. We can choose to center one of the mean vectors at zeroby subtracting the mean term ν̄ from each of the n elements of ν, which centers
those mean terms, and then adding this same term ν̄ to each of the p elements ofξ, which results in the same matrix of mean terms as the uncentered terms.
The last non-identifiability constraint we must deal with are the residual varianceterms, ψi,j = θi ηj . When p > n we look to see whether the largest element of θ
minus the smallest element of θ is greater than 3.0; if so, we scale θ by a constantso that θmax − θmin = 3.0, and we scale η by this same constant. Similarly, whenn > p, we check that ηmax − ηmin ≤ 3.0, and, if not, we rescale the η vector sothat this equation holds, and we scale θ by this same constant for consistency. We
recognize that this is not an ideal way to induce identifiability, but it appears togive reasonable solutions in practice.

5In some situations with certain data sets, we found SFA to converge to a spikein the likelihood, where a single individual’s genotype is modeled almost perfectlyby the factors (with a correspondingly low residual variance term). If this happens, there might be an issue with non-uniformity of variance across the loci (inwhich case, removing alleles with very low MAF may help). Also trying different
random starting points might find a good solution. Another possibility is to putthe variance terms on the loci instead of the individuals, or using SFAm. Finally,
if the same individual continues to be modeled repeatedly with different random
starting points, this person might be removed from the data set.
Because we choose not to update the expected values of the loading matrix Λ
between the CM steps, monotone convergence of the log marginal likelihood is notguaranteed, although in practice it appears to converge well. We find that convergence is reached for the applications described here after fewer than 500 iterations.
For each genotype data set, we run SFA multiple times with random seeds, setting the number of factors as described in the text; results presented in figures
are a representative example. A C++ package containing the general SFA code isavailable for download at http://stephenslab.uchicago.edu/software.html,
along with some data sets for application.