Restricted maximum-likelihood method for learning latent variance components in gene expression data with known and unknown confounders

Abstract Random effects models are popular statistical models for detecting and correcting spurious sample correlations due to hidden confounders in genome-wide gene expression data. In applications where some confounding factors are known, estimating simultaneously the contribution of known and latent variance components in random effects models is a challenge that has so far relied on numerical gradient-based optimizers to maximize the likelihood function. This is unsatisfactory because the resulting solution is poorly characterized and the efficiency of the method may be suboptimal. Here, we prove analytically that maximum-likelihood latent variables can always be chosen orthogonal to the known confounding factors, in other words, that maximum-likelihood latent variables explain sample covariances not already explained by known factors. Based on this result, we propose a restricted maximum-likelihood (REML) method that estimates the latent variables by maximizing the likelihood on the restricted subspace orthogonal to the known confounding factors and show that this reduces to probabilistic principal component analysis on that subspace. The method then estimates the variance–covariance parameters by maximizing the remaining terms in the likelihood function given the latent variables, using a newly derived analytic solution for this problem. Compared to gradient-based optimizers, our method attains greater or equal likelihood values, can be computed using standard matrix operations, results in latent factors that do not overlap with any known factors, and has a runtime reduced by several orders of magnitude. Hence, the REML method facilitates the application of random effects modeling strategies for learning latent variance components to much larger gene expression datasets than possible with current methods.


Supplementary Figures
Figure S1: A. Number of hidden covariates inferred by LVREML as a function of the parameter ρ (the targeted total amount of variance explained by the known and hidden covariates), with θ (the minimum variance explained by a known covariate) set to retain 0, 5, 10, or 20 known covariates (genotype PCs) in the model. B. Same as panel A, with θ set to retain 50, 100, 150, or 200 genotype PCs in the model. The saturation of the number of hidden covariates with decreasing ρ for models with 100, 150, and 200 known covariates is a visual indicator that some of the dimensions in the linear subspace spanned by the known covariates do not explain sufficient variation in the expression data, and the relevance or possible redundancy of (some of) the known covariates for explaining variation in the expression data needs to be reconsidered.

S1 Preliminary results
In the sections below, we will repeatedly use the following results. The first result concerns linear transformations of normally distributed variables and can be found in most textbooks on statistics or probability theory: Lemma 1. Let x ∈ R n be a random, normally distributed vector, with µ ∈ R n , and Ψ ∈ R n×n a positive definite covariance matrix. For any linear transformation y = Mx with M ∈ R m×n , we have If the linear transformation y = Mx in this Lemma is overdetermined, that is, if m > n, then the transformed covariance matrix Ψ ′ = MΨM T will have a lower rank n than its dimension m, that is, Ψ ′ ∈ R m×m is a positive semi-definite matrix (i.e., has one or more zero eigenvalues). Thus we can extend the definition of normal distributions to include degenerate distributions with positive semi-definite covariance matrix, by interpreting them as the distributions of overdetermined linear combinations of normally distributed vectors. A degenerate one-dimensional normal distribution is simply defined as a δ-distribution, that is, for which can be derived as a limit σ 2 → 0 of normal distribution density functions N (µ, σ 2 ).
The second result is one that is attributed to von Neumann [1]: Let P, Q ∈ R n×n be two positive definite matrices. Then where π 1 ≥ · · · ≥ π n and χ 1 ≥ · · · ≥ χ n are the ordered eigenvalues of P and Q, respectively, and equality in eq. (S1) is achieved if and only if the eigenvector of P corresponding to π i is equal to the eigenvector of Q corresponding to χ n−i+1 , i = 1, . . . , n.

S2 The model
We will use the following notation: • Y ∈ R n×m is a matrix of gene expression data for m genes in n samples. The ith column of Y is denoted y i ∈ R n and corresponds to the vector of expression values for gene i. We assume that the data in each sample are centred, ∑ m i=1 y i = 0 ∈ R n . • Z ∈ R n×d is a matrix of values for d known confounders in the same n samples. The kth column of Z is denoted z k ∈ R n and corresponds to the data for confounding factor k.
• X ∈ R n×p is a matrix of values for p latent variables to be determined in the same n samples. The jth column of X is denoted x j ∈ R n .
To identify the hidden correlation structure of the expression data, we assume a linear relationship between expression levels and the known and latent variables, with random noise added: where v i ∈ R d and w i ∈ R p are jointly normally distributed random vectors, is a positive semi-definite matrix; the errors ϵ i ∈ R n are assumed to be independent and normally distributed, Note that our aim is to identify variance components shared across genes, and hence σ 2 is assumed to be the same for all i. By assumption, the errors are also independent of the effect sizes, and hence we can write By Lemma 1, y i is normally distributed with distribution where and we used the notation ⟨u, v⟩ = u T v to denote the inner product between two vectors in R n .
Defining matrices V ∈ R d×m and W ∈ R p×m , whose columns are the random effect vectors v i and w i , respectively, eq. (S2) can be written in matrix notation as Under the assumption that the columns y i of Y are independent samples of the distribution (S5), the likelihood of observing Y given covariate data Z, (unknown) latent variable data X and values for the hyper-parameters Θ = {σ 2 , A, B, D}, is given by Note that in standard mixed-model calculations, the distribution (S5) is often arrived at by integrating out the random effects. This is equivalent to application of Lemma 1.
To conclude, the log-likelihood is, upto an additive constant, and divided by half the number of genes: is the empirical covariance matrix.

S3 Systematic effects on the mean
Eq. (S2) only considers random effects, which leads to a model for studying systematic effects on the covariance between samples. We could also include fixed effects to model systematic effects on mean expression level. However, by centering the data, ∑ m i=1 y i = 0, the maximumlikelihood estimate of such fixed effects is always zero. To see this, let T ∈ R n×c be a matrix of c covariates with fixed effects β ∈ R c shared across genes (we are only interested in discovering systematic biases in the data). Then the minus log-likelihood (2) becomes Optimizing with respect to β leads to the equation

S4 Solution of the model without latent variables
We start by considering the problem of finding the maximum-likelihood solution in the absence of any latent variables, i.e. minimizing eq. (2) with with respect to B and σ 2 .
Note first of all that we may assume the set of confounding factors {z 1 , . . . , z d } to be linearly independent, because if not, the expression in eq. (S2) can be rearranged in terms of a linearly independent subset of factors whose coefficients are still normally distributed due to elementary properties of the multivariate normal distribution, see for instance the proof of Lemma 5 below. Linear independence of {z 1 , . . . , z d } implies that we must have d ≤ n and rank(Z) = d.
The singular value decomposition allows to decompose Z as Z = UΓV T , where U ∈ R n×n , Note that unitarity of U implies U T 1 U 2 = 0. Denote by H Z the space spanned by the columns (i.e. covariate vectors) of Z. The projection matrix P Z onto H Z is given by Using the basis of column vectors of U, we can write any matrix M ∈ R n×n as a partitioned matrix The following results for partitioned matrices are derived easily or can be found in [2]: Using this notation, the following result solves the model without latent variables: where λ min (·) denotes the smallest eigenvalue of a matrix. Then the maximum-likelihood solution subject to B being positive semi-definite and σ 2 ≥ 0, is given bŷ Proof. Using eq. (S7), we can write Hence, in the block matrix notation (S8), we have It follows that and, using eqs. (S10) and (S11), Applying Lemma 2 to the term tr(K −1 11 C 11 ), it follows that for the minimizerK,K 11 must have eigenvalues κ 1 ≥ · · · ≥ κ d with the same eigevectors u 1 , . . . , u d as C 11 . Expressing the minus log-likelihood in terms of these eigenvalues results in Minimizing with respect to the parameters κ i and σ 2 (i.e., setting their derivatives to zero) results in the solutionκ i = λ i for all i andσ 2 = tr(C 22 ) n−d . In other words,K 11 has the same eigenvalues and eigenvectors as C 11 , that is, This equation is satisfied ifB Eq. (S12) is a condition on the amount of variation in Y explained by the confounders Z, with λ min (C 11 ) being (proportional to) the minimum amount of variation explained by any of the dimensions spanned by the columns of Z, and 1 n−d tr(C 22 ) being the average amount of variation explained by the dimensions orthogonal to the columns of Z. Failure of this condition simply means that there must be other, latent variables that explain more variation than the known ones, which is precisely what we are seeking to detect.
A useful special case of Theorem 1 occurs when the number of confounders equals one. In this case, we are seeking maximum-likelihood solutions for K of the form where z ∈ R n is the confounding data vector. Let γ 2 = ∥z∥ 2 and u = 1 γ z. Then P z = uu T is the projection matrix onto z, C 11 = ⟨u, Cu⟩, and tr(C 22 ) = tr((1 − P z )C) = tr(C) − ⟨u, Cu⟩. By Theorem 1, we haveβ

S5 Solution of the model without known covariates
Next, consider a model without known covariates, i.e. with posterior sample covariance matrix This model is equivalent to probabilistic principal component analysis [3,4], and its maximumlikelihood solution is given by the first p eigenvectors or principal components with largest eigenvalues of C. Here we present a more direct proof of this fact than what can be found in the literature.

Lemma 3.
Without loss of generality, we may assume that the latent variables have unit norm, are linearly independent, and are mutually orthogonal.
Proof. If the latent variables do not have unit norm, define Next assume that the latent variables are not linearly independent, i.e. that rank(K X ) = r < p. Because K X is a symmetric matrix, we must have K X = ∑ r l=1 t l t T l for some set of linearly independent vectors t l ∈ R n . Define α ′ l = ∥t l ∥ and x ′ l = t l /∥t l ∥. Then x ′ l has unit norm and . Because we may now assume that rank(X) = p, and because α j > 0 for all j, the matrix XA 1 2 has singular value decomposition We will also need the following simple result: Lemma 4. Let λ 1 ≥ λ 2 ≥ · · · ≥ λ n > 0 be a decreasing sequence of positive numbers, and let 1 ≤ p < n. If there exists j > p such that λ p > λ j , then (S17) Proof. Eq. (S17) follows from because each term on the r.h.s. is non-negative, and at least one is strictly positive. Theorem 2. Let C ∈ R n×n be a positive definite matrix with eigenvalues λ 1 ≥ · · · ≥ λ n and corresponding eigenvectors u 1 , . . . , u n , and let either p = n or 1 ≤ p < n such that there exists j > p with λ p > λ j . Then the maximum-likelihood solution Proof. By Lemma 3, we can assume that X has orthonormal columns, and hence there exist V ∈ R n×(n−p) such that Q = (X, V) ∈ R n×n is unitary, Q T Q = QQ T = 1. Hence K = XAX T + σ 2 1 has the spectral decomposition and hence where v l ∈ R n are the columns of V.
Assume that the α 2 j are ordered, α 2 1 ≥ · · · ≥ α 2 p . Applying von Neumann's Lemma 2 gives with equality if and only if x j = u j for j = 1, . . . , p v l = u p+l for l = 1, . . . , n − p Hence, independent of the values for α j , the maximum-likelihood latent variables are the eigenvectors of C corresponding to the p largest eigenvalues. Minimizing eq, (S18) w.r.t. α 2 j and σ 2 then gives By Lemma 4, α 2 j > 0 for all j.
Note that plugging the maximum-likelihood values in the likelihood function gives Either p can be set a priori small enough such that condition (S17) is satisfied, or else the value of p with smallest L min satisfying this condition can be found easily from eq. (S19).
Note also that in the models of [3,4], uniform prior variances are assumed (α 2 1 = · · · = α 2 p = 1), such that X is defined upto an arbitrary rotation, because XX T = (XR)(XR) T for any rotation matrix R. In our model, there is no such rotational freedom (if A is assumed to be diagonal), except if C has eigenvalues with multiplicities greater than one, when there is some freedom to choose the corresponding eigenvectors.

S6.1 Orthogonality of known and hidden confounders
Lemma 5. Without loss of generality, we may assume that the latent variables are orthogonal to the known confounders: Proof. As in Section S4, let P Z again be the projection matrix on the space spanned by the known covariates z k (i.e. the columns of Z). For any choice of latent variables x j , we have for some matrix of linear coefficients M = (m kq ) ∈ R d×p , and with ⟨s k ,x j ⟩ = 0 for all k. Or, in matrix notation Plugging this in eq. (S2), results in and hence, using Lemma 1, it follows that This is still of exactly the same form as eq. (S4). Hence model (S21) is identical to model (S2), but has hidden covariates orthogonal to the known covariates.
Note that we can parameterize the model with hidden variables orthogonal to the known confounders, Z T X = 0, but only if we allow the covariances of their effects on gene expression, Cov(v i , w i ) = D, to be non-zero. Equivalently, we can parameterize the model such that the random effects of hidden variables are statistically independent of the effects of the known confounders, Cov(v i , w i ) = 0, but only if we allow the hidden variables to overlap with the known confounders, Z T X ̸ = 0. Mathematically, the choice of orthogonal hidden factors will be much more convenient.
Note also that a transformation to orthogonal hidden factors always induces non-zero covariances among the known confounders via the term MAM T . Hence an important difficulty with the model where B is assumed to be diagonal, as used in [5], comes from the fact that non-orthogonal hidden variables are needed to model off-diagonal covariances between the known confounders. It is much more intuitive to model these directly by assuming a general covariance matrix.

S6.2 Restricted maximum-likelihood solution for the latent variables
Lemma 6. Without loss of generality, we may assume that the latent variables have unit norm, are linearly independent, and are mutually orthogonal.
Proof. The proof is identical to the proof of Lemma 3 -it is straightforward to verify that the transformation to orthonormal variables also do not change the form of the off-diagonal term ZDX T in the covariance matrix K, but merely lead to a reparameterization of the matrix D.
To solve the full model, we follow an approach similar to the standard restricted maximumlilelihood method for linear mixed models [6,7]: we write the negative log-likelihood function L = log det(K) + tr(K −1 C) as a sum where L 2 will be the log-likelihood restricted to the subspace orthogonal to the known confounders Z. We will estimate the latent variables X and their effect covariances A by maximizing L 2 , and estimate the effect covariances B and D involving the known confounders by maximizing L 1 . Solving for the latent variables on a restricted subspace is motivated by the observation that if y ∈ R n is a sample from the model (S2), that is, p(y) = N (0, K), then In other words, restricted to the subspace orthogonal to Z, the general model becomes a probablistic PCA model where all variation in the data is explained by the latent variables.
Using standard results for the marginal and conditional distributions of a multivariate Gaussian, we have , where we used the partitioned matrix notation of eq. (S8). In particular, , and hence L 2 depends only on X, A and σ 2 . The restricted maximum likelihood solution for the latent variables follows immediately: where the minimum is taken over all X with X T Z = 0, and all positive semi-definite diagonal matriceŝ A. If there exists j > p such that λ p > λ j , then where λ 1 ≥ λ 2 ≥ · · · ≥ λ n−d are the sorted eigenvalues of C 22 with corresponding eigenvectors w 1 , . . . , w n−d ∈ R n−d , and W p = (w 1 , . . . , w p ) ∈ R (n−d)×p is the matrix with the first p eigenvectors of C 22 as columns.
Proof. DefiningX = U T 2 X ∈ R (n−d)×p , we have K 22 =XAX T + σ 2 1, and L 2 becomes precisely the minus log-likelihood of the model without known covariates (Section S5), as a function of the latent variablesX on the reduced (n − d)-dimensional space orthogonal to the known confounders Z. Hence by Theorem 2, where λ 1 ≥ λ 2 ≥ · · · ≥ λ n−d are the sorted eigenvalues of C 22 and W p ∈ R (n−d)×p is the matrix having the corresponding first p eigenvectors as columns. Note thatÂ is positive semi-definite by Lemma 4 and the assumption that there exists j > p such that λ p > λ j . It remains to 'pull-back'X to the original n-dimensional space, using the orthogonality condition (S20):X This proves eqs. (S23) and (S24).

S6.3 Solution for the variance parameters given the latent variables
WithX,Â andσ 2 determined by the minimization of L 2 in Theorem 3, L 2 (X,Â,σ 2 ) is constant in terms of the parameters B and D that remain to be optimized. Hence optimizing L 1 with respect to these parameters is the same as optimizing the total negative log-likelihood L(X,Â, B, D,σ 2 ) w.r.t. B and D. We have: where as before is the singular value decomposition of Z, and W p = (w 1 , . . . , w p ) ∈ R (n−d)×p is the matrix with the first p eigenvectors of C 22 as columns.
Proof. Note that the conditions B and B − DÂ −1 D T positive semi-definite are to ensure that the matrix B D D TÂ is positive semi-definite. Next note that withX T known, the covariance matrix K can be written as Hence the total log-likelihood is identical to the model with known covariatesZ = ZX and no latent variables (Section S4). The unconstrained maximizing solution (that is, where A and σ 2 are also optimized) for the model with known covariatesZ is given by Theorem 1.
Due toX T Z = 0 and the definition ofX, the singular value decomposition ofZ is given bỹ where the columns of U 3 ∈ R n×(n−d−p) span the space orthogonal to the columns ofZ.
Working out the block matrix product results in: Hence, also the estimateÂ ′ =Â. Because the unconstrained optimization of L givenX results in the same estimate for A and σ 2 as the intial constrained optimization where these parameters were given, it follows that also the estimates of B and D must be the same:

S6.4 LVREML maximizes the variance explained
It is tempting to ask whether the combined solution from Theorems 3 and 4 optimizes the total likelihood among all possible p-dimensional sets of latent variables. To address this problem, let X ∈ R n×p be an arbitrary matrix of latent variables whose columns are normalized, mutually orthogonal and orthogonal to the columns of Z, X T X = 1 and X T Z = 0. Because U 2 is only defined upto a rotation, we can always choose Then Plugging these values into the negative log-likelihood function results in a function that depends only on X: Proposition 2. Let X ∈ R n×p be an arbitrary choice of latent variables with associated maximumlikelihood estimates for the covariance parameters given by Proposition 1. Then, upto an additive constant Proof. Recall from Theorem 2 that the maximum-likelihood estimate for K given X and its associated maximum-likelihood parameters estimates is given bŷ while the covariance matrix C can be written as Using equation (S11) for the determinant of a partitioned matrix, we have Ignoring the constants log det(U T 1 CU 1 ) and n which do not depend on X, we obtain eq. (S29).
Due to the determinant term in eq. (S29), it is not clear whether the restricted maximumlikelihood solutionX of Theorem 3 (with its associated maximum-likelihood covariance parameters of Theorem 4) is the absolute minimizer of L X , However, we do have the following result: Theorem 5. The restricted maximum-likelihood solutionX of Theorem 3 is the set of p latent variables that minimizes the residual variance among all choices of p latent variables, Proof. By Proposition 1 and the arguments leading up to it, we can write where as before C 22 = U T 2 CU 2 is the restriction of C to the (n − d)-dimensional subspace orthogonal to the d known covariates, and the columns of U T 2 X and U T 2 Q span mutually orthogonal subspaces within this (n − d)-dimensional space. Hence (n − d − p)σ 2 (X) = tr(Q T CQ T ) is the trace of C 22 over the residual (n − d − p)-dimensional space orthogonal to the latent variables, within the subspace orthogonal to the d known covariates. By the Courant-Fisher min-max theorem for eigenvalues [2], the (n − d − p)-dimensional subspace of R n−d with smallest trace is the subspace spanned by the eigenvectors of C 22 corresponding to its (n − d − p) smallest eigenvalues. By Theorem 3, this is exactly the subspace obtained by choosing X equal to the restricted maximum-likelihood solutionX.

S7 Selecting covariates and the latent dimension
Two practical problems remain: how to choose the latent variable dimension parameter p and which known covariates to include?
To choose p, we will use the following result: tr(C) = tr(K) = tr(ZBZ T ) + tr(XÂX T ) + nσ 2 Proof. Use Theorem 4 to compute where the last step uses the cyclical property of the trace and the fact that U T Hence tr(K) = tr(ZBZ) + tr(XÂX) + nσ 2 = tr(C 11 ) + tr(C 22 ) = tr(C) Because C = (YY T )/m, the eigenvalues of C are (proportional to) the squared singular values of the expression data Y. Hence tr(ZBZ)/ tr(C) is the proportion of variation in Y explained by the known covariates, tr(XÂX)/ tr(C) the proportion of variation explained by the latent variables, and nσ 2 / tr(C) is the residual variance.
Our method for determining the number of latent variables lets the user decide a priori the minimum amount of variation ρ in the data that should be explained by the known and latent confounders. It follows that given ρ, a "target" value for σ 2 is where the minimum with λ min (C 11 ) is taken to ensure that of condition (S26) remains valid. Because the eigenvalues λ 1 , . . . , λ n−d are sorted, the function increases with decreasing p. Hence given ρ, we definep aŝ that is, we choosep to be the smallest number of latent variables that explain at least a proporition of variation ρ of Y, while guaranteeing that the conditions for all mathematical results derived in this document are valid.
Note that unless all eigenvalues of C 22 are identical,p always exists. Once the desired number of latent variablesp is defined, the latent factorsX, the variance parametersÂ, and the residual variance estimateσ 2 (which will be the largest possible value less than or equal to the target value σ 2 (ρ)) are determined by Theorem 3. Once those are determined, the remaining covariance parametersB andD are determined by Theorem 4.
A second practical problem occurs when the rank of Z exceeds the number of samples, such that any subset of n linearly independent covariates explains all of the variation in Y. To select a more relevant subset of covariates, we rapidly screen all candidate covariates using the model with a single known covariate (Section S4) to compute the varianceβ 2 explained by that covariate alone (eq. (S16)). We then keep only those covariates for whichβ 2 ≥ θ tr(C), where θ > 0 is the second free parameter of the method, namely the minimum amount of variation explained by a known covariate on its own. The selected covariates are ranked according to their value ofβ 2 , and a linearly independent subset is generated, starting from the covariates with highestβ 2 .

S8 Downstream analyses
The inferred maximum-likelihood hidden factorsX and sample covariance matrixK are typically used to create a dataset of residuals corrected for spurious sample correlations, to increase the power for detecting eQTLs, or as data-derived endophenotypes [5,8]. We briefly review these tasks and how they compare between LVREML and PANAMA hidden factors.

S8.1 Correcting data for spurious sample correlations
To remove spurious correlations due to the known and latent variance components from the expression data Y ∈ R n×m (see Section S2), the residualsŷ i ∈ R n for gene i with original data y i (a column of Y) are contructed aŝ where the variance parameters σ 2 c,i and σ 2 e,i are fit separately for each gene i [5]. Hence two solutions for the latent factors that give rise to the sameK (as observed in Section 2.3 for LVREML and PANAMA) will result in the same residuals.

S8.2 Adjusting for known and latent covariates in eQTL association analyses
Two approaches for mapping eQTLs are commonly used in this context. The first approach tests for an association between SNP s j and gene y i using a mixed model, where the SNP is treated as a fixed effect, constructing likelihood ratio statistics as LOD i,j = log N (y i | θs j , σ 2 c,iK + σ 2 e,i 1) N (y i | 0, σ 2 c,iK + σ 2 e,i 1) , where the variance parameters σ 2 c,i and σ 2 e,i are fit separately for each gene i [5]. Hence for latent factor solutions that give rise to the sameK the association analyses will again be identical.
The second approach performs a linear regression of a gene's expression data, typically using the corrected dataŷ i , on the SNP genotypes s j , using the known and inferred factors as covariates [8], that is, a linear model is fit wherê where Z andX are the matrices of known and estimated latent factors, respectively, and a i ∈ R d and b i ∈ R p are their respective regression coefficients.
Since maximum-likelihood solutions for the hidden factors by LVREML and PANAMA differ by a linear combination with the known factors Z that transforms models with hidden factors orthogonal to Z to equivalent models with hidden factors overlapping with Z, and vice versa (see Section S6.1), it is clear that the same linear transformation will also result in equivalent linear association models in eq. (S30). Hence this type of analysis will also be equivalent between the hidden factors inferred by both approaches.

S8.3 Mapping the genetic architecture of latent variables
Inferred latent variables are sometimes treated as endophenotypes whose genetic architecture is of interest. In this case SNPs are identified that are strongly associated with the latent variables. Different solutions for the latent variables will then clearly result in different sets of significantly associated SNPs.
Using the maximum-likelihood LVREML inferred latent variables that are orthogonal to known confounders is advantageous in this context, because • The LVREML latent variables are uniquely defined. All other solutions that give rise to the same covariance matrix estimateK can be written as a linear combination of the known covariates and the LVREML covariates (see Section S6.1).
• When interpreting associated SNPs, there is no risk of attributing biological meaning to a latent variable that is due to the signal coming from the overlapping known covariates.
To remove the dependence of genetic association analyses on the choice of equivalent sets of latent variables, we recommend performing a multi-trait GWAS on the joint set of known and latent confounders. If the standard multivariate association test based on canonical correlation analysis [9] is used, results will again be identical between equivalent choices of latent variables, because together with the known confounders they all span the same linear subspace.