Unifying and Generalizing Methods for Removing Unwanted Variation Based on Negative Controls

TL;DR

This paper introduces RUV*, a unifying framework for removing unwanted variation in gene expression studies, which clarifies relationships between existing methods and enables new approaches like RUVB based on Bayesian analysis.

Contribution

The paper presents a general framework, RUV*, that unifies and extends existing methods for unwanted variation removal, facilitating development of new algorithms.

Findings

RUV* clarifies connections between existing RUV methods.

RUVB, based on RUV*, performs competitively in simulations.

Challenges remain in achieving consistent calibration across datasets.

Abstract

Unwanted variation, including hidden confounding, is a well-known problem in many fields, particularly large-scale gene expression studies. Recent proposals to use control genes --- genes assumed to be unassociated with the covariates of interest --- have led to new methods to deal with this problem. Going by the moniker Removing Unwanted Variation (RUV), there are many versions --- RUV1, RUV2, RUV4, RUVinv, RUVrinv, RUVfun. In this paper, we introduce a general framework, RUV*, that both unites and generalizes these approaches. This unifying framework helps clarify connections between existing methods. In particular we provide conditions under which RUV2 and RUV4 are equivalent. The RUV* framework also preserves an advantage of RUV approaches --- their modularity --- which facilitates the development of novel methods based on existing matrix imputation algorithms. We illustrate this by implementing RUVB, a version of RUV* based on Bayesian factor analysis. In realistic simulations based on real data we found that RUVB is competitive with existing methods in terms of both power and calibration, although we also highlight the challenges of providing consistently reliable calibration among data sets.

Equations

$${\mathbf{Y}}_{n\times p}={\mathbf{X}}_{n\times k}{\mathbf{\beta }}_{k\times p}+{\mathbf{Z}}_{n\times q}{\mathbf{\alpha }}_{q\times p}+{\mathbf{E}}_{n\times p},$$

$${\mathbf{Q}}^{⊺}\mathbf{Y}=\mathbf{R}\mathbf{\beta }+{\mathbf{Q}}^{⊺}\mathbf{Z}\mathbf{\alpha }+{\mathbf{Q}}^{⊺}\mathbf{E}.$$

$$\displaystyle\boldsymbol{R}_{1}=\left(\begin{array}[]{cc}\boldsymbol{R}_{11}&\boldsymbol{R}_{12}\\ 0&\boldsymbol{R}_{22}\end{array}\right).$$

$${\mathbf{Y}}_1$$

$${\mathbf{Y}}_2$$

$${\mathbf{Y}}_3$$

$${\beta }_{ij}=0\text{ for all }i=k_1+1,\dots ,k,\text{ and }j\in \mathcal{C},$$

$${\mathbf{Y}}_{2\mathcal{C}}$$

$$={\mathbf{Z}}_2{\hat{\mathbf{\alpha }}}_{\mathcal{C}}+{\mathbf{E}}_{2\mathcal{C}},$$

$$e_{2\mathcal{C}ij}$$

$${\hat{\mathbf{Z}}}_2={\mathbf{Y}}_{2\mathcal{C}}{\hat{\mathbf{\alpha }}}_{\mathcal{C}}^{⊺}({\hat{\mathbf{\alpha }}}_{\mathcal{C}}{\hat{\mathbf{\alpha }}}_{\mathcal{C}}^{⊺}{)}^{-1}.$$

$${\hat{\mathbf{Z}}}_2={\mathbf{Y}}_{2\mathcal{C}}{\hat{\mathbf{\Sigma }}}_{\mathcal{C}}^{-1}{\hat{\mathbf{\alpha }}}_{\mathcal{C}}^{⊺}({\hat{\mathbf{\alpha }}}_{\mathcal{C}}{\hat{\mathbf{\Sigma }}}_{\mathcal{C}}^{-1}{\hat{\mathbf{\alpha }}}_{\mathcal{C}}^{⊺}{)}^{-1}.$$

$${\hat{\mathbf{\beta }}}_2={\mathbf{R}}_{22}^{-1}({\mathbf{Y}}_2-{\hat{\mathbf{Z}}}_2\hat{\mathbf{\alpha }}).$$

$$\hat{\mathbf{\beta }}=({\mathbf{X}}^{⊺}\mathbf{S}\mathbf{X}{)}^{-1}{\mathbf{X}}^{⊺}\mathbf{S}\mathbf{Y},$$

$$\hat{\mathbf{\alpha }}$$

$$\hat{\mathbf{\Sigma }}$$

$${\hat{\mathbf{\beta }}}_2={\mathbf{R}}_{22}^{-1}({\mathbf{Y}}_2-{\hat{\mathbf{Z}}}_2\hat{\mathbf{\alpha }}).$$

$${\mathcal{F}}_1(\mathbf{Y})$$

$${\mathcal{F}}_2(\mathbf{Y})$$

$${\text{RUV2}}_{old}({\mathcal{F}}_2)={\text{RUV2}}_{new}({\mathcal{F}}_1).$$

$$\displaystyle\boldsymbol{W}\left(\begin{array}[]{c}\boldsymbol{Q}_{2}^{\intercal}\\ \boldsymbol{Q}_{3}^{\intercal}\end{array}\right)\boldsymbol{Y}=\boldsymbol{W}\left(\begin{array}[]{c}\boldsymbol{Q}_{2}^{\intercal}\\ \boldsymbol{Q}_{3}^{\intercal}\end{array}\right)\boldsymbol{X}_{2}\boldsymbol{\beta}_{2}+\boldsymbol{W}\left(\begin{array}[]{c}\boldsymbol{Q}_{2}^{\intercal}\\ \boldsymbol{Q}_{3}^{\intercal}\end{array}\right)\boldsymbol{Z}\boldsymbol{\alpha}+\boldsymbol{W}\left(\begin{array}[]{c}\boldsymbol{Q}_{2}^{\intercal}\\ \boldsymbol{Q}_{3}^{\intercal}\end{array}\right)\boldsymbol{E}$$

$$\displaystyle\boldsymbol{W}\left(\begin{array}[]{c}\boldsymbol{Q}_{2}^{\intercal}\\ \boldsymbol{Q}_{3}^{\intercal}\end{array}\right)\boldsymbol{Y}=\boldsymbol{W}\left(\begin{array}[]{c}\boldsymbol{R}_{22}\\ \boldsymbol{0}\end{array}\right)\boldsymbol{\beta}_{2}+\boldsymbol{W}\left(\begin{array}[]{c}\boldsymbol{Q}_{2}^{\intercal}\\ \boldsymbol{Q}_{3}^{\intercal}\end{array}\right)\boldsymbol{Z}\boldsymbol{\alpha}+\boldsymbol{W}\left(\begin{array}[]{c}\boldsymbol{Q}_{2}^{\intercal}\\ \boldsymbol{Q}_{3}^{\intercal}\end{array}\right)\boldsymbol{E}.$$

$$\displaystyle\tilde{\boldsymbol{Y}}\mathrel{\mathop{\mathchar 58\relax}}=\boldsymbol{W}\left(\begin{array}[]{c}\boldsymbol{Q}_{2}^{\intercal}\\ \boldsymbol{Q}_{3}^{\intercal}\end{array}\right)\boldsymbol{Y},\ \tilde{\boldsymbol{X}}\mathrel{\mathop{\mathchar 58\relax}}=\boldsymbol{W}\left(\begin{array}[]{c}\boldsymbol{R}_{22}\\ \boldsymbol{0}\end{array}\right),\ \tilde{\boldsymbol{Z}}\mathrel{\mathop{\mathchar 58\relax}}=\boldsymbol{W}\left(\begin{array}[]{c}\boldsymbol{Q}_{2}^{\intercal}\\ \boldsymbol{Q}_{3}^{\intercal}\end{array}\right)\boldsymbol{Z},\text{ and }\tilde{\boldsymbol{E}}\mathrel{\mathop{\mathchar 58\relax}}=\boldsymbol{W}\left(\begin{array}[]{c}\boldsymbol{Q}_{2}^{\intercal}\\ \boldsymbol{Q}_{3}^{\intercal}\end{array}\right)\boldsymbol{E}.$$

$$\tilde{\mathbf{Y}}=\tilde{\mathbf{X}}\mathbf{\beta }+\tilde{\mathbf{Z}}\mathbf{\alpha }+\tilde{\mathbf{E}}.$$

$$\{\hat{\mathbf{\Sigma }}({\mathbf{Q}}^{⊺}\mathbf{Y}),\hat{\mathbf{Z}}({\mathbf{Q}}^{⊺}\mathbf{Y})\mathbf{A}(\mathbf{Y}),\mathbf{A}(\mathbf{Y}{)}^{-1}\hat{\mathbf{\alpha }}({\mathbf{Q}}^{⊺}\mathbf{Y})\}=\{\hat{\mathbf{\Sigma }}(\mathbf{Y}),{\mathbf{Q}}^{⊺}\hat{\mathbf{Z}}(\mathbf{Y}),\hat{\mathbf{\alpha }}(\mathbf{Y})\},$$

$${\text{RUV2}}_{old}(\mathcal{F})={\text{RUV2}}_{new}(\mathcal{F}).$$

$${\mathcal{F}}_2(\mathbf{Y})$$

$$=\{\hat{\mathbf{\Sigma }}(\mathbf{Y}),\mathbf{Q}{\mathbf{Q}}^{⊺}\hat{\mathbf{Z}}(\mathbf{Y}),\hat{\mathbf{\alpha }}(\mathbf{Y})\}$$

$$=\{\hat{\mathbf{\Sigma }}(\mathbf{Y}),\hat{\mathbf{Z}}(\mathbf{Y}),\hat{\mathbf{\alpha }}(\mathbf{Y})\}={\mathcal{F}}_1(\mathbf{Y}).$$

$${\text{RUV2}}_{old}(\mathcal{F})={\text{RUV2}}_{old}({\mathcal{F}}_2)={\text{RUV2}}_{new}({\mathcal{F}}_1)={\text{RUV2}}_{new}(\mathcal{F}).$$

$$\hat{\mathbf{Z}}(\mathbf{Y})$$

$$\hat{\mathbf{\alpha }}(\mathbf{Y})$$

$$\hat{\mathbf{\Sigma }}(\mathbf{Y})$$

$${\mathbf{Y}}_{2\mathcal{C}}$$

$${\mathbf{Y}}_{2\bar{\mathcal{C}}}$$

$${\mathbf{Y}}_{3\mathcal{C}}$$

Full text

††footnotetext: MSC 2010 subject classifications: primary 62J15, 62F15; secondary 62H25, 62P10

Keywords and phrases: unobserved confounding, hidden confounding, correlated tests, negative controls, RNA-seq, gene expression, unwanted variation, batch effects

1 Introduction

Many experiments and observational studies in genetics are overwhelmed with unwanted sources of variation. Examples include: processing date [Akey et al., 2007], the lab that collected a sample [Irizarry et al., 2005], the batch in which a sample was processed [Leek et al., 2010], and subject attributes such as environmental factors [Gibson, 2008] and ancestry [Price et al., 2006]. These factors, if ignored, can result in disastrously wrong conclusions [Gilad and Mizrahi-Man, 2015]. They can induce dependencies between samples, and inflate test statistics, making it difficult to control false discovery rates [Efron, 2004, 2008, 2010].

Many of the sources of variation mentioned above are likely to be observed, in which case standard methods exist to control for them [Johnson et al., 2007]. However, every study likely also contains unobserved sources of unwanted variation, and these can cause equally profound problems [Leek and Storey, 2007] — even in the ideal case of a randomized experiment. To illustrate this we took $20$ samples from an RNA-seq dataset [GTEx Consortium, 2015] and randomly assigned them into two groups of 10 samples. Since group assignment is entirely independent of the expression levels of each gene, the group labels are theoretically unassociated with all genes and any observed “signal” must be artefactual. Figure 1 shows histograms of the $p$-values from two-sample $t$-tests for three different randomizations. In each case the distribution of the $p$-values differs greatly from the theoretical uniform distribution. Thus, even in this ideal scenario where group labels were randomly assigned, problems can arise. One way to understand this is to note that the same randomization is being applied to all genes. Consequently, if many genes are affected by an unobserved factor, and this factor happens by chance to be correlated with the randomization, then the $p$-value distributions will be non-uniform. In this sense the problems here can be viewed as being due to correlation among the $p$ values; see Efron [2010] for extensive discussion. (The issue of whether the problems in any given study are caused by correlation, confounding, or something different is both interesting and subtle; see discussion in Efron [2010], Schwartzman [2010] for example. For this reason we adopt the “unwanted variation” terminology from Gagnon-Bartsch and Speed [2012], rather than alternative terminologies such as “hidden confounding”.)

In recent years many methods have been introduced to try to solve problems due to unwanted variation. Perhaps the simplest approach is to estimate sources of unwanted variation using principal components analysis [Price et al., 2006], and then to control for these factors by using them as covariates in subsequent analyses. Indeed, in genome-wide association studies this simple method is widely used. However, in gene expression studies it suffers from the problem that the principal components will typically also contain the signal of interest, so controlling for them risks removing that signal. To address this Leek and Storey [2007, 2008] introduced Surrogate Variable Analysis (SVA), which uses an iterative algorithm to attempt to estimate latent factors that do not include the signal of interest (see also Lucas et al. [2006] for related work). To account for unwanted variation, SVA assumes a factor-augmented regression model (see Section 2.1), from which there is an old literature [Fisher and Mackenzie, 1923, Cochran, 1943, Williams, 1952, Tukey, 1962, Gollob, 1968, Mandel, 1969, 1971, Efron and Morris, 1972, Freeman et al., 1973, Gabriel, 1978, and others]. Since SVA, a large number of different approaches have emerged along similar lines, including Behzadi et al. [2007], Carvalho et al. [2008], Kang et al. [2008a, b], Stegle et al. [2008], Friguet et al. [2009], Kang et al. [2010], Listgarten et al. [2010], Stegle et al. [2010], Wu and Aryee [2010], Gagnon-Bartsch and Speed [2012], Fusi et al. [2012], Stegle et al. [2012], Sun et al. [2012], Gagnon-Bartsch et al. [2013], Mostafavi et al. [2013], Perry and Pillai [2013], Yang et al. [2013], Wang et al. [2015], Chen and Zhou [2016], among others.

As noted above, a key difficulty in adjusting for unwanted variation in expression studies is distinguishing between the effect of a treatment and the effect of factors that are correlated with a treatment. Available methods deal with this problem in different ways. In this paper we focus on the subset of methods that use “negative controls” to help achieve this goal. In the context of a gene expression study, a negative control is a gene whose expression is assumed a priori to be unassociated with all covariates (and treatments) of interest. Under this assumption, negative controls can be used to separate sources of unwanted variation from the treatment effects. The idea of using negative controls in this way appears in Lucas et al. [2006], and has been recently popularized by Gagnon-Bartsch and Speed [2012], Gagnon-Bartsch et al. [2013] in a series of methods and software going by the moniker Removing Unwanted Variation (RUV). There are several different methods, including RUV2 (for RUV 2-step), RUV4, RUVinv (a special case of RUV4), RUVrinv, RUVfun, and RUV1.

Understanding the relative merits and properties of the different RUV methods, which are all aimed at solving essentially the same problem, is a non-trivial task. The main contribution of this paper is to outline a general framework, RUV*, that encompasses all versions of RUV (Section 5). RUV* represents the problem as a general matrix imputation procedure, both providing a unifying conceptual framework, and opening up new approaches based on the large literature in matrix imputation. Our RUV* framework also provides a simple and modular way to account for uncertainty in the estimated sources of unwanted variation, which is an issue ignored by most methods. On the way to this general framework we make detailed connections between RUV2 and RUV4, exploiting the formulation in Wang et al. [2015].

The structure of the paper is as follows. We first review the formulation of Wang et al. [2015] (Section 2.1) in the context of RUV4 (Section 2.2). We then extend this to include RUV2 (Section 3), and develop necessary and sufficient conditions for a procedure to be a version of both RUV2 and RUV4 (Section 4). We call the resulting procedure RUV3. We then develop the general RUV* framework (Section 5), and illustrate it by implementing an approach based on Bayesian factor analysis that we call RUVB. In Section 6 we briefly discuss the issue of variance estimation, which turns out to be an important issue that can greatly affect empirical performance. In Section 7 we use realistic simulations, based on real data, to compare many of the available methods that use negative controls. We show that RUVB has competitive power and calibration compared with other methods, particularly when there are few control genes. We finish with a discussion in Section 8.

On notation: throughout we denote matrices using bold capital letters ($\boldsymbol{A}$), except for $\boldsymbol{\alpha}$ and $\boldsymbol{\beta}$, which are also matrices. Bold lowercase letters are vectors ($\boldsymbol{a}$), and non-bold lowercase letters are scalars ($a$). Where there is no chance for confusion, we use non-bold lowercase to denote scalar elements of vectors or matrices. For example, $a_{ij}$ is the $(i,j)$th element of $\boldsymbol{A}$ and $a_{i}$ is the $i$th element of $\boldsymbol{a}$. The notation $\boldsymbol{A}_{n\times m}$ denotes that the matrix $\boldsymbol{A}$ is an $n$ by $m$ matrix. The matrix transpose is denoted $\boldsymbol{A}^{\intercal}$ and the matrix inverse is denoted $\boldsymbol{A}^{-1}$. Sets are generally denoted with calligraphic letters ($\mathcal{A}$), and the complement of a set is denoted with a bar ($\bar{\mathcal{A}}$). Finally, $p(\boldsymbol{a}|\boldsymbol{b})$ denotes a probability density function of $\boldsymbol{a}$ conditional on $\boldsymbol{b}$.

2 RUV4

2.1 Review of the Two-step Rotation Method

Most existing approaches to this problem [Leek and Storey, 2007, 2008, Gagnon-Bartsch and Speed, 2012, Sun et al., 2012, Gagnon-Bartsch et al., 2013, Wang et al., 2015] use a low-rank matrix to capture unwanted variation. Specifically, they assume:

[TABLE]

where, in the context of a gene-expression study, $y_{ij}$ is the normalized expression level of the $j$th gene on the $i$th sample, $\boldsymbol{X}$ contains the observed covariates, $\boldsymbol{\beta}$ contains the coefficients of $\boldsymbol{X}$, $\boldsymbol{Z}$ is a matrix of unobserved factors (sources of unwanted variation), $\boldsymbol{\alpha}$ contains the coefficients of $\boldsymbol{Z}$, and $\boldsymbol{E}$ contains independent (Gaussian) errors with means 0 and column-specific variances $\operatorname*{var}(e_{ij})=\sigma_{j}^{2}$. In this model, the only known quantities are $\boldsymbol{Y}$ and $\boldsymbol{X}$.

To fit (1), it is common to apply a two-step approach (e.g. Gagnon-Bartsch et al. [2013], Sun et al. [2012], Wang et al. [2015]). The first step regresses out $\boldsymbol{X}$ and then, using the residuals of this regression, estimates $\boldsymbol{\alpha}$ and the $\sigma_{j}$’s. The second step then assumes that $\boldsymbol{\alpha}$ and the $\sigma_{j}$’s are known and estimates $\boldsymbol{\beta}$ and $\boldsymbol{Z}$. Wang et al. [2015] helpfully frame this two-step approach as a rotation followed by estimation in two independent models. We now review this approach.

First, we let $\boldsymbol{X}=\boldsymbol{Q}\boldsymbol{R}$ denote the QR decomposition of $\boldsymbol{X}$, where $\boldsymbol{Q}\in\mathbb{R}^{n\times n}$ is an orthogonal matrix ($\boldsymbol{Q}^{\intercal}\boldsymbol{Q}=\boldsymbol{Q}\boldsymbol{Q}^{\intercal}=\boldsymbol{I}_{n}$) and $\boldsymbol{R}_{n\times k}={{\boldsymbol{R}_{1}}\choose{0}}$, where $\boldsymbol{R}_{1}\in\mathbb{R}^{k\times k}$ is an upper-triangular matrix. Multiplying (1) on the left by $\boldsymbol{Q}^{\intercal}$ yields

[TABLE]

Suppose that $k=k_{1}+k_{2}$, where the first $k_{1}$ covariates of $\boldsymbol{X}$ are not of direct interest, but are included because of various modeling decisions (e.g. an intercept term, or covariates that need to be controlled for). The last $k_{2}$ columns of $\boldsymbol{X}$ are the variables of interest whose putative associations with $\boldsymbol{Y}$ the researcher wishes to test. Let $\boldsymbol{Y}_{1}\in\mathbb{R}^{k_{1}\times p}$ be the first $k_{1}$ rows of $\boldsymbol{Q}^{\intercal}\boldsymbol{Y}$, $\boldsymbol{Y}_{2}\in\mathbb{R}^{k_{2}\times p}$ be the next $k_{2}$ rows of $\boldsymbol{Q}^{\intercal}\boldsymbol{Y}$, and $\boldsymbol{Y}_{3}\in\mathbb{R}^{(n-k)\times p}$ be the last $n-k$ rows of $\boldsymbol{Q}^{\intercal}\boldsymbol{Y}$. Conformably partition $\boldsymbol{Q}^{\intercal}\boldsymbol{Z}$ into $\boldsymbol{Z}_{1}$, $\boldsymbol{Z}_{2}$, and $\boldsymbol{Z}_{3}$, and $\boldsymbol{Q}^{\intercal}\boldsymbol{E}$ into $\boldsymbol{E}_{1}$, $\boldsymbol{E}_{2}$, and $\boldsymbol{E}_{3}$. Let

[TABLE]

Finally, partition $\boldsymbol{\beta}={{\boldsymbol{\beta}_{1}}\choose{\boldsymbol{\beta}_{2}}}$ so that $\boldsymbol{\beta}_{1}\in\mathbb{R}^{k_{1}\times p}$ contains the coefficients for the first $k_{1}$ covariates and $\boldsymbol{\beta}_{2}\in\mathbb{R}^{k_{2}\times p}$ contains the coefficients for the last $k_{2}$ covariates. Then (2) may be written as three models

[TABLE]

Importantly, the error terms in (6), (7), and (8) are mutually independent. This follows from the easily-proved fact that $\boldsymbol{E}$ is equal in distribution to $\boldsymbol{Q}^{\intercal}\boldsymbol{E}$. The two-step estimation procedure mentioned above becomes: first, estimate $\boldsymbol{\alpha}$ and the $\sigma_{j}$’s using (8); second, estimate $\boldsymbol{\beta}_{2}$ and $\boldsymbol{Z}_{2}$ given $\boldsymbol{\alpha}$ and the $\sigma_{j}$’s using (7). Equation (6) contains the nuisance parameters $\boldsymbol{\beta}_{1}$ and is ignored.

2.2 Review of RUV4

One approach to distinguishing between unwanted variation and effects of interest is to use “control genes” [Lucas et al., 2006, Gagnon-Bartsch and Speed, 2012]. A control gene is a gene that is assumed a priori to be unassociated with the covariate(s) of interest. More formally, the set of control genes, $\mathcal{C}\subseteq\{1,\ldots,p\}$, has the property that

[TABLE]

and is a subset of the truly null genes. Examples of control genes used in practice are spike-in controls [Jiang et al., 2011] used to adjust for technical factors (such as sample batch) and housekeeping genes [Eisenberg and Levanon, 2013] used to adjust for both technical and biological factors (such as subject ancestry).

RUV4 [Gagnon-Bartsch et al., 2013] uses control genes to estimate $\boldsymbol{\beta}_{2}$ in the presence of unwanted variation. Let $\boldsymbol{Y}_{2\mathcal{C}}\in\mathbb{R}^{k_{2}\times m}$ denote the submatrix of $\boldsymbol{Y}_{2}$ with columns that correspond to the $m$ control genes. Similarly subset the relevant columns to obtain $\boldsymbol{\beta}_{2\mathcal{C}}\in\mathbb{R}^{k_{2}\times m}$, $\boldsymbol{\alpha}_{\mathcal{C}}\in\mathbb{R}^{q\times m}$, and $\boldsymbol{E}_{2\mathcal{C}}\in\mathbb{R}^{k_{2}\times m}$. The steps for RUV4, including a variation from Wang et al. [2015], are presented in Procedure 1. (For simplicity we focus on point estimates of effects here, deferring assessment of standard errors to Section 6.)

The key idea in Procedure 1 is that for the control genes model (7) becomes

[TABLE]

The equality in (10) follows from the property of control genes that $\boldsymbol{\beta}_{2\mathcal{C}}=\boldsymbol{0}$. Step 2 of Procedure 1 uses (10) to estimate $\boldsymbol{Z}_{2}$.

Step 1 of Procedure 1 requires a factor analysis of $\boldsymbol{Y}_{3}$. We formally define a factor analysis as follows.

Definition 1**.**

A factor analysis, $\mathcal{F}$, of rank $q\leq\min(n,p)$ on $\boldsymbol{Y}\in\mathbb{R}^{n\times p}$ is a set of three functions $\mathcal{F}=\{\hat{\boldsymbol{\Sigma}}(\boldsymbol{Y}),\hat{\boldsymbol{Z}}(\boldsymbol{Y}),\hat{\boldsymbol{\alpha}}(\boldsymbol{Y})\}$ such that $\hat{\boldsymbol{\Sigma}}(\boldsymbol{Y})\in\mathbb{R}^{p\times p}$ is diagonal with positive diagonal entries, $\hat{\boldsymbol{Z}}(\boldsymbol{Y})\in\mathbb{R}^{n\times q}$ has rank $q$, and $\hat{\boldsymbol{\alpha}}(\boldsymbol{Y})\in\mathbb{R}^{q\times p}$ has rank $q$.

RUV4 allows the analyst to use any factor analysis they desire. Thus, RUV4 is not a single method, but a collection of methods indexed by the factor analysis used. When we want to be explicit about this indexing, we will write RUV4$(\mathcal{F})$.

3 RUV2

RUV2 is a method for removing unwanted variation presented in Gagnon-Bartsch and Speed [2012] (predating RUV4 in Gagnon-Bartsch et al. [2013]). One contribution of our paper is to present RUV2 within the same rotation framework as RUV4 (Section 2.1).

Procedure 2 summarizes RUV2 as presented in Gagnon-Bartsch and Speed [2012]. The method is simple: first estimate the factors causing unwanted variation from the control genes, and then include these factors as covariates in the regression models for the non-control genes. However, this procedure does not deal with nuisance parameters. To deal with these, Gagnon-Bartsch et al. [2013] introduce an extension, also called RUV2. This extension first rotates $\boldsymbol{Y}$ and $\boldsymbol{X}$ onto the orthogonal complement of the columns of $\boldsymbol{X}$ corresponding to the nuisance parameters [equation (64) in Gagnon-Bartsch et al., 2013] and then applies Procedure 2.

Like RUV4, RUV2 is actually a class of methods indexed by the factor analysis used (Definition 1). We denote the class of RUV2 methods presented in Gagnon-Bartsch and Speed [2012] by $\text{RUV2}_{old}(\mathcal{F})$.

In Procedure 3 we present a class of methods, denoted $\text{RUV2}_{new}(\mathcal{F})$, that we will prove to be equivalent to $\text{RUV2}_{old}$.

Theorem 1**.**

For a given orthogonal matrix $\boldsymbol{Q}\in\mathbb{R}^{n\times n}$ and an arbitrary non-singular matrix $\boldsymbol{A}(\boldsymbol{Y})$ that (possibly) depends on $\boldsymbol{Y}$, suppose

[TABLE]

Then

[TABLE]

That is, Procedure 2 using factor analysis (19) is equivalent to Procedure 3 using factor analysis (20).

Proof.

For simplicity, we first assume that there are no nuisance covariates. Then ${\boldsymbol{Y}_{2\mathcal{C}}\choose\boldsymbol{Y}_{3\mathcal{C}}}=\boldsymbol{Q}^{\intercal}\boldsymbol{Y}$, where $\boldsymbol{Q}$ is the orthogonal matrix from the QR decomposition of $\boldsymbol{X}$ (Section 2.1). Thus, ${\hat{\boldsymbol{Z}}_{2}\choose\hat{\boldsymbol{Z}}_{3}}$ from Procedure 3 results from applying $\mathcal{F}_{1}$ on $\boldsymbol{Q}^{\intercal}\boldsymbol{Y}$ while $\hat{\boldsymbol{Z}}$ from Procedure 2 results from applying $\mathcal{F}_{2}$ on $\boldsymbol{Y}$. From the definitions of (19) and (20), we thus have that $\hat{\boldsymbol{Z}}$ from step 1 of Procedure 2 is in the same column space as $\boldsymbol{Q}{\hat{\boldsymbol{Z}}_{2}\choose\hat{\boldsymbol{Z}}_{3}}$ from step 1 of Procedure 3. $\hat{\boldsymbol{\beta}}_{2}$ from (18) contains the partial regression coefficients of $\boldsymbol{X}$ when including $\boldsymbol{Q}{\hat{\boldsymbol{Z}}_{2}\choose\hat{\boldsymbol{Z}}_{3}}$ as nuisance covariates (to show this, just calculate the MLE’s of $\boldsymbol{\beta}_{2}$ and $\boldsymbol{\alpha}$ using (7) and (8)). The estimates of $\boldsymbol{\beta}_{2}$ in step 2 of Procedure 2 are also partial regression coefficients of $\boldsymbol{X}$ when including $\hat{\boldsymbol{Z}}$ as nuisance covariates. Since the partial regression coefficients in Procedure 2 are only a function of $\hat{\boldsymbol{Z}}$ through its column space, and the partial regression coefficients in Procedure 3 are only a function of $\boldsymbol{Q}{\hat{\boldsymbol{Z}}_{2}\choose\hat{\boldsymbol{Z}}_{3}}$ through its column space, and these column spaces are the same, we have completed the proof for the case of no nuisance parameters.

To deal with nuisance parameters, Gagnon-Bartsch et al. [2013] rotate $\boldsymbol{X}$ and $\boldsymbol{Y}$ onto the orthogonal complement of the columns of $\boldsymbol{X}$ corresponding to the nuisance parameters prior to applying Procedure 2. If we partition $\boldsymbol{Q}=(\boldsymbol{Q}_{1},\boldsymbol{Q}_{2},\boldsymbol{Q}_{3})$ and $\boldsymbol{X}=(\boldsymbol{X}_{1},\boldsymbol{X}_{2})$, this is equivalent to using the model

[TABLE]

where $\boldsymbol{W}$ is some arbitrary (but known) $n-k_{1}$ by $n-k_{1}$ orthogonal matrix. Or, using the QR decomposition of $\boldsymbol{X}$, (30) is equal to

[TABLE]

Let

[TABLE]

Then, (30) is equal to

[TABLE]

We now just apply the arguments of the previous paragraph to (49), where there are no nuisance parameters and where the QR decomposition of $\tilde{\boldsymbol{X}}$ is just $\tilde{\boldsymbol{X}}=\boldsymbol{W}{\boldsymbol{R}_{22}\choose\boldsymbol{0}}$. This completes the proof. ∎

The equivalence of $\text{RUV2}_{\text{old}}$ and $\text{RUV2}_{\text{new}}$ in Theorem 1 involves using different factor analyses in each procedure. One can ask under what conditions the two procedures would be equivalent if given the same factor analysis. We now prove that it suffices for the factor analysis to be left orthogonally equivariant.

Definition 2**.**

A factor analysis of rank $q$ on $\boldsymbol{Y}\in\mathbb{R}^{n\times p}$ is left orthogonally equivariant if

[TABLE]

for all fixed orthogonal $\boldsymbol{Q}\in\mathbb{R}^{n\times n}$ and an arbitrary non-singular $\boldsymbol{A}(\boldsymbol{Y})\in\mathbb{R}^{q\times q}$ that (possibly) depends on $\boldsymbol{Y}$.

Corollary 1**.**

Suppose $\mathcal{F}$ is a left orthogonally equivariant factor analysis. Then

[TABLE]

Proof.

Suppose $\mathcal{F}_{2}$ from (20) is left orthogonally equivariant, then

[TABLE]

Let $\mathcal{F}\mathrel{\mathop{\mathchar 58\relax}}=\mathcal{F}_{1}=\mathcal{F}_{2}$. From the results of Theorem 1, we have that

[TABLE]

∎

A well-known example of a factor analysis that is left orthogonally equivariant is a truncated singular value decomposition (SVD). That is, let $\boldsymbol{Y}=\boldsymbol{U}\boldsymbol{D}\boldsymbol{V}^{\intercal}$ be the SVD of $\boldsymbol{Y}$. Let $\boldsymbol{D}^{(q)}=\operatorname*{diag}(d_{11},\ldots,d_{qq},0,\ldots,0)\in\mathbb{R}^{n\times n}$ be a diagonal matrix whose first $q$ diagonal elements are the same as in $\boldsymbol{D}$ and the last $n-q$ diagonal elements are [math]. Then

[TABLE]

for any $\pi\in[0,1]$. (Without loss of generality, for the remainder of this paper, we let $\pi=1$.) Notably, this factor analysis is the only option in the R package ruv [Gagnon-Bartsch, 2015].

For the rest of this paper, we use RUV2 to refer to Procedure 3 and not Procedure 2, even if the factor analysis used is not orthogonally equivariant. (By Theorem 1, this corresponds to Procedure 2 with some other factor analysis.)

4 RUV3

Gagnon-Bartsch et al. [2013] provide a lengthy heuristic discussion comparing RUV2 with RUV4 (their section 3.4). However, they provide no mathematical equivalencies. In this section, we introduce RUV3, a procedure that is both a version of RUV2 and a version RUV4. We show that it is the only such procedure that is both RUV2 and RUV4. RUV3 uses a partitioned matrix imputation procedure from Owen and Wang [2016] to estimate the hidden factors. The coefficients are then estimated as in RUV2 and RUV4.

4.1 The RUV3 Procedure

The main goal in all methods is to estimate $\boldsymbol{\beta}_{2\bar{\mathcal{C}}}$, the coefficients corresponding to the non-control genes. This involves incorporating information from four models

[TABLE]

We can rearrange the rows and columns of ${\boldsymbol{Y}_{2}\choose\boldsymbol{Y}_{3}}$ to write these models in matrix form:

[TABLE]

The major difference between RUV2 and RUV4 is how the estimation procedures interact in (67); see Figure 2 for a pictorial representation. RUV2 performs factor analysis on ${\boldsymbol{Y}_{2\mathcal{C}}\choose\boldsymbol{Y}_{3\mathcal{C}}}$, then regresses $\boldsymbol{Y}_{3\bar{\mathcal{C}}}$ on the estimated factor loadings. RUV4 performs factor analysis on $(\boldsymbol{Y}_{3\mathcal{C}},\boldsymbol{Y}_{3\bar{\mathcal{C}}})$, then regresses $\boldsymbol{Y}_{2\mathcal{C}}$ on the estimated factors. The main goal in both, however, is to estimate $\boldsymbol{Z}_{2}\boldsymbol{\alpha}_{\bar{\mathcal{C}}}$ given $\boldsymbol{Y}_{2\mathcal{C}}$, $\boldsymbol{Y}_{3\mathcal{C}}$, and $\boldsymbol{Y}_{3\bar{\mathcal{C}}}$. It is not clear why one should prefer either the rows or the columns to perform a factor analysis first, then perform a respective regression on the columns or rows.

We can frame the estimation of $\boldsymbol{Z}_{2}\boldsymbol{\alpha}_{\bar{\mathcal{C}}}$ as a matrix imputation problem where the missing values are a submatrix. That is, we try to estimate $\boldsymbol{Z}_{2}\boldsymbol{\alpha}_{\bar{\mathcal{C}}}$ given only the values of $\boldsymbol{Y}_{2\mathcal{C}}$, $\boldsymbol{Y}_{3\mathcal{C}}$, and $\boldsymbol{Y}_{3\bar{\mathcal{C}}}$. In the context of matrix imputation (and not removing unwanted variation), Owen and Wang [2016], generalizing the methods of Owen and Perry [2009] to the case of heteroscedastic noise, suggest that after applying a factor analysis to $\boldsymbol{Y}_{3\mathcal{C}}$, one use the estimates

[TABLE]

and then set $\widehat{\boldsymbol{Z}_{2}\boldsymbol{\alpha}_{\bar{\mathcal{C}}}}=\hat{\boldsymbol{Z}}_{2}\hat{\boldsymbol{\alpha}}_{\bar{\mathcal{C}}}$. This corresponds to a factor analysis followed by two regressions followed by an imputation step. Following the theme of this paper, we would add an additional step and estimate $\boldsymbol{\beta}_{2\bar{\mathcal{C}}}$ with

[TABLE]

In this section, we prove that this estimation procedure, presented in Procedure 4, is a version of both RUV2 (Section 4.3) and RUV4 (Section 4.2). Indeed, it is the only such procedure that is both a version of RUV2 and RUV4 (Section 4.4). As such, we call it RUV3 (also presented pictorially in Figure 2). Like RUV2 and RUV4, RUV3 is a class of methods indexed by the factor analysis used. We sometimes explicitly denote this indexing by $\text{RUV3}(\mathcal{F})$.

4.2 Connection to RUV4

The astute reader will have noticed that (71) is the same as (13). RUV3 can be viewed as a version of RUV4 with a particular factor analysis. Specifically, any factor analysis applied during RUV4 of the following form will result in RUV3:

[TABLE]

Or, more simply, RUV4 is equal to RUV3 if, in RUV4 one uses any factor analysis such that $\hat{\boldsymbol{\alpha}}_{\bar{\mathcal{C}}}=(\hat{\boldsymbol{Z}}_{3}^{\intercal}\hat{\boldsymbol{Z}}_{3})^{-1}\hat{\boldsymbol{Z}}_{3}^{\intercal}\boldsymbol{Y}_{3\bar{\mathcal{C}}}$ and $\hat{\boldsymbol{\Sigma}}_{\mathcal{C}}$, $\hat{\boldsymbol{Z}}_{3}$, and $\hat{\boldsymbol{\alpha}}_{\mathcal{C}}$ are functions of $\boldsymbol{Y}_{3\mathcal{C}}$ but not $\boldsymbol{Y}_{3\bar{\mathcal{C}}}$

Actually, using a truncated SVD on $(\boldsymbol{Y}_{3\mathcal{C}},\boldsymbol{Y}_{3\bar{\mathcal{C}}})$ (equations (56) to (58)) results in the equalities (76) to (78) if one ignores the functional dependencies. That is, using the truncated SVD on $(\boldsymbol{Y}_{3\mathcal{C}},\boldsymbol{Y}_{3\bar{\mathcal{C}}})$ one can easily prove the relationships

[TABLE]

The main difference, then, between RUV3 and RUV4 as implemented in the ruv R package [Gagnon-Bartsch, 2015] is that in RUV3 $\hat{\boldsymbol{Z}}_{3}$ has a functional dependence only on $\boldsymbol{Y}_{3\mathcal{C}}$ and not $\boldsymbol{Y}_{3\bar{\mathcal{C}}}$.

4.3 Connection to RUV2

Similar to the relationship between RUV3 and RUV4, RUV3 may also be viewed as a version of RUV2 with a particular factor analysis. Specifically any factor analysis applied during RUV2 of the following form will result in RUV3:

[TABLE]

In simpler terms, RUV2 is the same as RUV3 if, in RUV2 one uses any factor analysis such that $\hat{\boldsymbol{Z}}_{2}=\boldsymbol{Y}_{2\mathcal{C}}\hat{\boldsymbol{\Sigma}}_{\mathcal{C}}^{-1}\hat{\boldsymbol{\alpha}}_{\mathcal{C}}^{\intercal}(\hat{\boldsymbol{\alpha}}_{\mathcal{C}}\hat{\boldsymbol{\Sigma}}_{\mathcal{C}}^{-1}\hat{\boldsymbol{\alpha}}_{\mathcal{C}}^{\intercal})^{-1}$ and $\hat{\boldsymbol{\alpha}}_{\mathcal{C}}$, $\hat{\boldsymbol{Z}}_{3}$, and $\boldsymbol{\boldsymbol{\Sigma}}_{\mathcal{C}}$ are functions of $\boldsymbol{Y}_{3\mathcal{C}}$ but not $\boldsymbol{Y}_{2\mathcal{C}}$.

Similar to Section 4.2, using the truncated SVD on ${\boldsymbol{Y}_{2\mathcal{C}}\choose\boldsymbol{Y}_{3\mathcal{C}}}$ — assuming homoscedastic variances rather than heteroscedastic variances — will result in equalities (82) to (86) except for the functional dependencies. That is, when using the truncated SVD on ${\boldsymbol{Y}_{2\mathcal{C}}\choose\boldsymbol{Y}_{3\mathcal{C}}}$ one can show that the following relationships hold:

[TABLE]

The main difference, then, between RUV2 and RUV3 is that in RUV3 $\hat{\boldsymbol{\alpha}}_{\mathcal{C}}$ has a functional dependence only on $\boldsymbol{Y}_{3\mathcal{C}}$ and not $\boldsymbol{Y}_{2\mathcal{C}}$.

One apparent disadvantage to using the original RUV2 pipeline (Procedure 2 or 3), is that after the factor analysis of step 1, one already has estimates of $\boldsymbol{\alpha}_{\mathcal{C}}$ and $\boldsymbol{\Sigma}_{\mathcal{C}}$, yet one ignores these estimates and re-estimates them in step 2 with $\hat{\boldsymbol{\alpha}}_{\mathcal{C}}=(\hat{\boldsymbol{Z}}_{3}^{\intercal}\hat{\boldsymbol{Z}}_{3})^{-1}\hat{\boldsymbol{Z}}_{3}^{\intercal}\hat{\boldsymbol{Y}}_{3C}$ and $\hat{\boldsymbol{\Sigma}}_{\mathcal{C}}=\operatorname*{diag}[(\boldsymbol{Y}_{3\mathcal{C}}-\hat{\boldsymbol{Z}}_{3}\hat{\boldsymbol{\alpha}}_{\mathcal{C}})^{\intercal}(\boldsymbol{Y}_{3\mathcal{C}}-\hat{\boldsymbol{Z}}_{3}\hat{\boldsymbol{\alpha}}_{\mathcal{C}})]/(n-q)$. In Procedures 2 or 3, the estimates from step 1 are in general different from the estimates from step 2 and it is not at all clear which estimates one should prefer. RUV3 obviates this problem by the constrained factor analysis (86).

4.4 RUV2 and RUV4 if and only if RUV3

We have shown in Sections 4.2 and 4.3 that RUV3 can be considered a variant of both RUV2 and RUV4. But the converse is easily proved to also be true.

Theorem 2**.**

Suppose a procedure is both a version of RUV4 (Procedure 1) and RUV2 (Procedure 3). Then it is also a version of RUV3 (Procedure 4).

Proof.

If the procedure is a version of RUV2, then (16) holds. But if the procedure is a version of RUV4, then (13) holds. These are two properties of RUV3 (equations (71) and (72)).

It remains to show that $\hat{\boldsymbol{Z}}_{3}$ and $\hat{\boldsymbol{\alpha}}_{\mathcal{C}}$ are functions only of $\boldsymbol{Y}_{3\mathcal{C}}$. But this is clear since if the procedure is RUV4, these quantities are functions only of $\boldsymbol{Y}_{3\mathcal{C}}$ and $\boldsymbol{Y}_{3\bar{\mathcal{C}}}$, while if the procedure is RUV2 these quantities are functions only of $\boldsymbol{Y}_{2\mathcal{C}}$ and $\boldsymbol{Y}_{3\mathcal{C}}$. Since the procedure is both RUV2 and RUV4, this necessarily implies that these quantities are functions only of $\boldsymbol{Y}_{3\mathcal{C}}$. ∎

To summarize, RUV3 can be viewed as both a version of RUV2 and a version of RUV4 and if a procedure is both a version of RUV2 and a version of RUV4 then it is a version of RUV3.

5 A more general framework: RUV*

A key insight that arises from our unification of RUV2 and RUV4 (and RUV3) into a single framework is that they share a common goal: estimation of $\boldsymbol{Z}_{2}\boldsymbol{\alpha}_{\bar{\mathcal{C}}}$, which represents the combined effects of all sources of unwanted variation on $\boldsymbol{Y}_{2\bar{\mathcal{C}}}$.

This insight suggests a more general approach to the problem: any matrix imputation procedure could be used to estimate $\boldsymbol{Z}_{2}\boldsymbol{\alpha}_{\bar{\mathcal{C}}}$ — RUV2, RUV3, and RUV4 are just three versions that rely heavily on linear associations between submatrices. Indeed, we need not even assume a factor model for the form of the unwanted variation. And we can further incorporate uncertainty in the estimates. In this section we develop these ideas to provide a more general framework for removing unwanted variation, which we call RUV*.

5.1 More general approaches to matrix imputation

To generalize RUV methods to allow for more general approaches to matrix imputation we generalize (67) to

[TABLE]

where $\boldsymbol{\Omega}$ is the unwanted variation parameterized by some $\boldsymbol{\phi}$. When the unwanted variation is represented by a factor model, we have that $\boldsymbol{\phi}=\{\boldsymbol{Z},\boldsymbol{\alpha}\}$ and $\boldsymbol{\Omega}(\boldsymbol{\phi})=\boldsymbol{Z}\boldsymbol{\alpha}$.

The simplest version of RUV* fits this model in two steps:

Use any appropriate procedure to estimate $\boldsymbol{\Omega}_{2\bar{\mathcal{C}}}(\boldsymbol{\phi})$ given $\{\boldsymbol{Y}_{2\mathcal{C}},\boldsymbol{Y}_{3\mathcal{C}},\boldsymbol{Y}_{3\bar{\mathcal{C}}}\}$; 2. 2.

Estimate $\boldsymbol{\beta}_{2}$ by

[TABLE]

This idea, and its relationship with other RUV approaches, are illustrated in Figure 3. Rather than restrict factors to be estimated via linear regression, RUV* allows any imputation procedure to be used to estimate $\boldsymbol{\Omega}_{2\bar{\mathcal{C}}}(\boldsymbol{\phi})$. This opens up a large literature on matrix imputation for use in removing unwanted variation with control genes [Hoff, 2007, Allen and Tibshirani, 2010, Candes and Plan, 2010, Stekhoven and Bühlmann, 2012, van Buuren, 2012, Josse et al., 2016, for example].

5.2 Relationship to RUVfun

Gagnon-Bartsch et al. [2013] describe a general framework they call RUVfun, for RUV-functional. In our notation, and within the rotated model framework of Section 2.1, RUVfun may be described as

Perform factor analysis on $(\boldsymbol{Y}_{3\mathcal{C}},\boldsymbol{Y}_{3\bar{\mathcal{C}}})$ to obtain an estimate $\hat{\boldsymbol{\alpha}}$, 2. 2.

Let $\hat{\boldsymbol{\alpha}}_{j}$ denote the $j$th column of $(\hat{\boldsymbol{\alpha}}_{\mathcal{C}},\hat{\boldsymbol{\alpha}}_{\bar{\mathcal{C}}})$ and $\boldsymbol{y}_{j}$ denote the $j$th column of $(\boldsymbol{Y}_{2\mathcal{C}},\boldsymbol{Y}_{2\bar{\mathcal{C}}})$. Train a function $f$ using responses $\boldsymbol{y}_{j}$ and predictors $\hat{\boldsymbol{\alpha}}_{j}$ for $j=1,\ldots,m$ (recall, we have $m$ control genes). That is, fit

[TABLE]

and call the resulting predictor $\hat{f}$. 3. 3.

Estimate $\boldsymbol{Y}_{2\bar{\mathcal{C}}}$ using predictors $\hat{\boldsymbol{\alpha}}_{\bar{\mathcal{C}}}$. That is,

[TABLE] 4. 4.

Estimate $\boldsymbol{\beta}_{2\bar{\mathcal{C}}}$ with $\boldsymbol{R}_{22}^{-1}(\boldsymbol{Y}_{2\bar{\mathcal{C}}}-\hat{\boldsymbol{Y}}_{2\bar{\mathcal{C}}})$.

This is the way it is presented in Section 3.8.2 of Gagnon-Bartsch et al. [2013], but typically they take the factor analysis step to mean just setting $\hat{\boldsymbol{\alpha}}=(\boldsymbol{Y}_{3\mathcal{C}},\boldsymbol{Y}_{3\bar{\mathcal{C}}})$. A pictorial representation of RUVfun is presented in the third panel of Figure 3. The only difference between the RUVfun diagram in Figure 3 and the RUV4 diagram in Figure 2 is that “regression” was changed to “train any function”.

RUVfun, though more general than RUV4, is a special case of RUV*. And RUV* is more general: for example, RUV2 is a version of RUV* but not of RUVfun. Indeed, RUV* generalizes RUVfun in three key ways. First, RUVfun uses only one column of $\hat{\boldsymbol{\alpha}}$ to estimate one column of $\boldsymbol{Y}_{2\bar{\mathcal{C}}}$ while RUV* allows for joint estimation of $\boldsymbol{Y}_{2\bar{\mathcal{C}}}$. Second, RUVfun assumes that each column of the rotated $\boldsymbol{Y}$ matrix is independent and identically distributed [Gagnon-Bartsch et al., 2013, page 41] while RUV* does not. Indeed, some matrix imputation approaches use column covariances to great effect [Allen and Tibshirani, 2010]. Third, RUVfun uses only $\boldsymbol{Y}_{3\mathcal{C}}$ to train the prediction function, whereas RUV* can use all elements in the rotated $\boldsymbol{Y}$ matrix.

5.3 Incorporating uncertainty in estimated unwanted variation

Like previous RUV methods, the second step of RUV* described above treats the estimate of $\boldsymbol{\Omega}_{2\bar{\mathcal{C}}}(\boldsymbol{\phi})$ from the first step as if it were “known” without error. Here we generalize this, using Bayesian ideas to propagate the uncertainty from the first step to the next.

Although the use of Bayesian methods in this context is not new [Stegle et al., 2008, 2010, Fusi et al., 2012, Stegle et al., 2012], our development here shares one of the great advantages of the RUV methods: modularity. That is, RUV methods separate the analysis into smaller self-contained steps: the factor analysis step and the regression step. Modularity is widely used in many fields: mathematicians modularize results using theorems, lemmas and corollaries; computer scientists modularize code using functions and classes. Modularity has many benefits, including: (1) it is easier to conceptualize an approach if it is broken into small simple steps, (2) it is easier to discover and correct mistakes, and (3) it is easier to improve an approach by improving specific steps. These advantages also apply to statistical analysis and methods development. For example, in RUV if one wishes to use a new method for factor analysis then this does not require a whole new approach — one simply replaces the truncated SVD with the new factor analysis.

To describe this generalized RUV* we introduce a latent variable $\tilde{\boldsymbol{Y}}_{2\bar{\mathcal{C}}}$ and write (97) as

[TABLE]

Now consider the following two-step procedure:

Use any appropriate procedure to obtain a conditional distribution $h(\tilde{\boldsymbol{Y}}_{2\bar{\mathcal{C}}})\mathrel{\mathop{\mathchar 58\relax}}=p(\tilde{\boldsymbol{Y}}_{2\bar{\mathcal{C}}}|\mathcal{Y}_{m})$, where $\mathcal{Y}_{m}\mathrel{\mathop{\mathchar 58\relax}}=\{\boldsymbol{Y}_{2\mathcal{C}},\boldsymbol{Y}_{3\mathcal{C}},\boldsymbol{Y}_{3\bar{\mathcal{C}}}\}$. 2. 2.

Perform inference for $\boldsymbol{\beta}_{2}$ using the likelihood

[TABLE]

where $\delta(\cdot)$ indicates the Dirac delta function.

Of course, in step 2 one could do classical inference for $\boldsymbol{\beta}_{2}$, or place a prior on $\boldsymbol{\beta}_{2}$ and perform Bayesian inference.

This procedure requires an analytic form for the conditional distribution $h$. An alternative is to assume that we can sample from this conditional distribution, which yields a convenient sample-based (or “multiple imputation”) RUV* algorithm.

Use any appropriate procedure to obtain samples $\tilde{\boldsymbol{Y}}_{2\bar{\mathcal{C}}}^{(1)},\ldots,\tilde{\boldsymbol{Y}}_{2\bar{\mathcal{C}}}^{(t)}$ from a conditional distribution $p(\tilde{\boldsymbol{Y}}_{2\bar{\mathcal{C}}}|\mathcal{Y}_{m})$. 2. 2.

Approximate the likelihood for $L(\boldsymbol{\beta}_{2})$ by using the fact that $\hat{\boldsymbol{\beta}}_{2}^{(i)}\mathrel{\mathop{\mathchar 58\relax}}=\boldsymbol{R}_{22}^{-1}(\boldsymbol{Y}_{2\bar{\mathcal{C}}}-\tilde{\boldsymbol{Y}}_{2\bar{\mathcal{C}}}^{(i)})$ are sampled from a distribution proportional to $L(\boldsymbol{\beta}_{2})$ [which is guaranteed to be proper; Appendix A.4].

For example, here in step 2 we can approximate the likelihood for each element of $\boldsymbol{\beta}_{2}$ by a normal likelihood

[TABLE]

where $\hat{\beta}_{2j}$ and $\hat{s}_{j}$ are respectively the mean and standard deviation of $\hat{\boldsymbol{\beta}}_{2}^{(i)}$. Alternatively, a $t$ likelihood can be used. Either approach provides an estimate and standard error for each element of ${\boldsymbol{\beta}}_{2}$ that accounts for uncertainty in the estimated unwanted variation. This is in contrast to the various methods used by the other RUV approaches to estimate the standard errors which do not account for this uncertainty (Section 6). Here we use these values to rank the “significance” of genes by the value of $\hat{\beta}_{2j}/\hat{s}_{j}$. They could also be used as inputs to the empirical Bayes method in Stephens [2016] to obtain measurements of significance related to false discovery rates.

Other approaches to inference in Step 2 are also possible. For example, given a specific prior on $\boldsymbol{\beta}_{2}$, Bayesian inference for $\boldsymbol{\beta}_{2}$ could be performed by re-weighting these samples according to this prior distribution (Appendix A.1). This re-weighting yields an arbitrarily accurate approximation to the posterior distribution $p(\boldsymbol{\beta}_{2}|\mathcal{Y}_{m},\boldsymbol{Y}_{2\bar{\mathcal{C}}})$ (Appendix A.2). Posterior summaries using this re-weighting scheme are easy to derive (Appendix A.5).

5.4 RUVB: Bayesian factor analysis in RUV*

To illustrate the potential for RUV* to produce new methods for removing unwanted variation we implemented a version of RUV*, using a Markov chain Monte Carlo scheme to fit a simple Bayesian Factor analysis model, and hence perform the sampling-based imputation in Step 1 of RUV*. See Appendix A.3 for details. We refer to this method as RUVB.

6 Estimating Standard Errors

For simplicity we have focused our descriptions of RUV2, RUV3, RUV4, and RUV* on point estimation for $\boldsymbol{\beta}_{2}$. In practice, to be useful, all of these methods must also provide standard errors for these estimates. Several different approaches to this problem exist, and we have found in empirical comparisons (e.g. Section 7 below) that the approach taken can greatly affect results, particularly calibration of interval estimates. In this section we therefore briefly review some of these approaches.

The simplest approach is to treat the estimated $\hat{\boldsymbol{Z}}$ as the true value of ${\boldsymbol{Z}}$, and then use standard theory from linear models to estimate the standard errors of the estimated coefficients. That is, first estimate $\hat{\boldsymbol{Z}}$ using any of the RUV approaches, then regress $\boldsymbol{Y}$ on $(\boldsymbol{X},\hat{\boldsymbol{Z}})$ and obtain estimated standard errors (for coefficients of $\boldsymbol{X}$) in the usual way. This is the default option in the ruv R package. The cate R package implements this (with the nc.var.correction and calibrate parameters both set to FALSE), but without the usual degrees of freedom correction in estimated standard errors. Though asymptotically this will not matter, we have found that for small sample sizes this can substantially hurt performance due to downward-biased standard errors.

Gagnon-Bartsch et al. [2013] noted that the standard errors estimated using this simple approach can be poorly behaved, essentially because the $\hat{\boldsymbol{Z}}$ are estimated and not known. They suggested several approaches to calibrating these standard errors using control genes. The approach that we found to work best in our comparisons (at least, when there are many control genes — see Section 7) is to multiply the estimated standard errors by a factor $\lambda$ (i.e. set $\tilde{s}_{i}=\lambda\hat{s}_{i}$) which is estimated from control genes by

[TABLE]

where $\hat{\beta_{j}}$ and $\hat{s_{j}}$ are the estimated coefficients and their standard errors (Equation (236) in Gagnon-Bartsch et al. [2013]). In our empirical comparisons below we refer to this procedure as “control gene calibration”. Gagnon-Bartsch et al. [2013] use heuristic arguments to motivate (111) in the context of studies with just one covariate of interest. In Appendix A.7, we extend (111) to the case when there is more than one covariate of interest and formally justify it with maximum likelihood arguments.

Sun et al. [2012] take a different approach to calibration, which does not use control genes, but is motivated by the assumption that most genes are null. Specifically they suggest centering the $t$-statistics $\hat{\beta}_{j}/\hat{s}_{j}$ by their median and scaling them by their median absolute deviation (MAD):

[TABLE]

The motivation for this adjustment is that if most genes are null, then normalizing by robust estimates of the null $t$-statistics’ center and scale will make the null $t$-statistics more closely match their theoretical distribution. This adjustment of $t$ statistics is closely connected with the variance calibration (111). Indeed, if we assume that the median of the $t$-statistics is approximately zero, then it is effectively equivalent to variance calibration with

[TABLE]

in place of (111).

In addition to MAD calibration, Wang et al. [2015] offer asymptotic arguments for an additive variance adjustment. This additive inflation term is particularly important when there are few control genes, and is specific to the RUV4 estimator (unlike (111) and (113) which can be applied to any estimator).

Finally, before development of RUV-like methods, the benefits of using empirical Bayes variance moderation (EBVM) [Smyth, 2004] were widely recognized in gene expression analyses. Variance moderation can be applied in combination with the variance calibration methods discussed above: for example, the ruv::variance_adjust function in R applies EBVM before applying (111). EBVM can similarly be incorporated into other methods. For RUV3 and CATE we can use EBVM either before or after the generalized least squares (GLS) step of their respective algorithms (equation (71) for RUV3 and equation (13) for CATE).

7 Empirical Evaluations

7.1 Simulation Approach

We now use simulations based on real data to compare methods (focusing on methods that use control genes, although the same simulations could be useful more generally). In brief, we use random subsets of real expression data to create “null data” that contains real (but unknown) “unwanted variation”, and then modify these null data to add (known) signal. We compare methods in their ability to reliably detect real signals and avoid spurious signals. Because they are based on real data, our simulations involve realistic levels of unwanted variation. However, they also represent a “best-case” scenario in which treatment labels were randomized with respect to the factors causing this unwanted variation. They also represent a best case scenario in that the control genes given to each method are simulated to be genuinely null. Even in this best-case scenario unwanted variation is a major issue, and, as we shall see, obtaining well calibrated inferences is challenging.

In more detail: we took the top 1000 expressed genes from the RNA-seq data on muscle samples from the Genotype Tissue Expression Consortium [GTEx Consortium, 2015]. For each dataset in our simulation study, we randomly selected $n$ samples ($n=6$, $10$, $20$, or $40$). We then randomly assigned half of these samples to be in one group and the other half to be in a second group. So our design matrix $\boldsymbol{X}\in\mathbb{R}^{n\times 2}$ contains two columns — a column of ones and a column that is an indicator for group assignment.

At this point, all genes are theoretically null, as in the datasets of our introduction (Section 1). We used this “all-null” scenario as one setting in our simulation studies (similar to the simulations in Rocke et al. [2015]). For other settings, we added signal to a randomly selected proportion of genes $\pi_{1}=1-\pi_{0}$ ($\pi_{1}=0.1$ or $0.5$). To add signal, we first sampled the effect sizes from a $N(0,0.8^{2})$. The standard deviation, $0.8$, was chosen by trial and error so that the classification problem would be neither too easy nor too hard. Let

[TABLE]

be the effect sizes, where $j_{\ell}\in\Omega$, the set of non-null genes. Then we drew a $\boldsymbol{W}$ matrix of the same dimension as our RNA-seq count matrix $\boldsymbol{Z}$ by

[TABLE]

where $\delta(z_{ij_{\ell}})$ indicates a point mass at $z_{ij_{\ell}}$. We then used $\boldsymbol{W}$ as our new RNA-seq matrix. Prior to running each method, we took a simple $\log_{2}$ transformation of the elements in $\boldsymbol{W}$ to obtain our $\boldsymbol{Y}$ matrix.

To motivate this approach, suppose $z_{ij}\sim\operatorname*{Poisson}(\lambda_{j})$, then

[TABLE]

Hence,

[TABLE]

So $a_{j}$ is approximately the $\log_{2}$-fold difference between the two groups.

7.2 Summary of Methods Compared

We compared standard ordinary least squares regression (OLS) against five other approaches: RUV2, RUV3, RUV4, CATE (the GLS variant of RUV4), and RUVB. We tried each of these effect-estimation methods with different approaches to variance estimation (Section 6). Specifically, for the non-RUVB methods we considered:

• •

Variance moderation (EBVM) versus no moderation.

• •

Variance calibration (both MAD and control-gene based) vs no calibration.

• •

For RUV3 and CATE: EBVM before GLS or after GLS.

• •

For CATE: additive variance inflation vs no additive variance inflation.

Altogether this gave 6 different OLS approaches, 6 RUV2 approaches, 9 RUV3 approaches, 6 RUV4 approaches, and 15 CATE approaches. (We did not implement CATE with additive variance inflation and EBVM before GLS because this implementation is not straightforward given the current cate software.)

For RUVB, we considered four approaches to producing mean and variance estimates:

• •

Using sample-based posterior summaries (Appendix A.5).

• •

Using the normal approximation to the likelihood in Equation (110).

• •

Using a $t$ approximation to the likelihood, replacing Equation (110) with a $t$ density with $n-k-q$ degrees of freedom.

• •

Using a $t$ approximation as above, followed by EBVM.

Additional technical considerations are discussed in Appendix A.6.

Although the large number of methods considered may seem initially excessive, we have found that there is often more variation in performance among different versions of a method than among the different families of method (RUV2, RUV3, RUV4/CATE, and RUVB). That is, choices among strategies such as EBVM and variance calibration may matter as much (or more) in practice as choices between families such as RUV2 vs RUV4.

7.3 Comparisons: Sensitivity vs Specificity

To compare methods in their ability to distinguish null and non-null genes we calculated the area under the receiver operating characteristic curve (AUC) for each method as the significance threshold is varied. (Multiplicative variance calibration does not affect AUC because it does not change the significance rankings among genes; thus we do not discuss variance calibration methods in this section.)

The clearest result here is that all the methods consistently outperform standard OLS (Supplementary Figure 1). This emphasizes the benefits of removing unwanted variation in improving power to detect real effects. For small sample size comparisons (e.g. 3 vs 3) the gains are smaller — though still apparent — presumably because reliably estimating the unwanted variation is harder for small samples.

A second clear pattern is that the use of variance moderation (EBVM) consistently improved AUC performance: the best-performing method in each family used EBVM. As might be expected, these benefits of EBVM are greatest for smaller sample sizes (Supplementary Figure 1).

Compared with these two clear patterns, differences among the best-performing methods in each family are more subtle. Figure 5 compares the AUC of the best method in each family with that of RUVB, which performed best overall in this comparison. (Results are shown for $\pi_{0}=0.5$; results for $\pi_{0}=0.9$ are similar). We highlight four main results:

RUVB has the best mean AUC among all methods we explored; 2. 2.

RUV4/CATE methods perform less well (relative to RUVB) when there are few control genes and the sample size is large; 3. 3.

In contrast, RUV2 methods perform less well (relative to RUVB) when the sample size is small and there are few control genes; 4. 4.

RUV3 performs somewhat stably (relative to RUVB) across the sample sizes.

7.4 Comparisons: Calibration

We also assessed the calibration of methods by examining the empirical coverage of their nominal $95\%$ confidence intervals for each effect (based on standard theory for the relevant $t$ distribution in each case). Because the variance calibration methods can have a strong effect (on all methods) we initially consider methods without variance calibration.

We begin by examining “typical” coverage for each method in each scenario by computing the median (across datasets) of the empirical coverage. We found that, without variance calibration, all method families except RUV4/CATE were able to achieve satisfactory typical coverage — somewhere between 0.94 and 0.97 — across all scenarios (Figure 6a) shows results for $\pi_{0}=0.5$; other values yielded similar results, not shown). The best performing method (in terms of typical coverage) in the RUV4/CATE family was CATE with only additive variance inflation. However, this method was often overly conservative in scenarios with few control genes, especially with larger sample sizes.

Although these median coverage results are encouraging, in practice having small variation in coverage among datasets is also important. That is, we would like methods to have near-95% coverage in most data-sets, and not only on average. Here the results (Figure 6b; Supplementary Figure 2) are less encouraging: the coverage of the methods with good typical coverage (median coverage close to 95%) varied considerably among datasets. This said, variability does improve for larger sample sizes and more control genes, and in this case all methods improve noticeably on OLS (Figure 6b, right facet). A particular concern is that, across all these methods, for many datasets, empirical coverage can be much lower than the nominal goal of 95%. Such datasets might be expected to lead to problems with over-identification of significant null genes (“false positives”), and under-estimation of false discovery rates when using either frequentist or Empirical Bayes FDR-related methodology [e.g. Benjamini and Hochberg, 1995, Storey, 2003, Stephens, 2016].

To summarize variability in coverage — as well as any tendency to be conservative or anti-conservative — we calculated the proportion of datasets where the actual coverage deviated substantially from $95\%$, which we defined as being either less than $90\%$ or more than $97.5\%$. Figure 7 shows these proportions for each method. Here we have also included methods that use variance calibration, as the results help highlight the effects of these calibration methods. The key findings are:

RUVB (the normal and sample-based versions) has “balanced” errors in coverage: its empirical coverage is as likely to be too high as too low. 2. 2.

MAD calibration tends to produce highly conservative coverage — that is, its coverage is very often much larger than the claimed 95%, and seldom much lower. This will tend to reduce false positive significant results, but also substantially reduce power to detect real effects. The exception is that when all genes are null ($\pi_{0}=1$), MAD calibration works well for larger sample sizes. These results are likely explained partly by non-null genes biasing upwards the variance calibration parameter, an issue also noted in Sun et al. [2012]. 3. 3.

Control-gene calibration is often anti-conservative when there are few control genes. However, it can work well when the sample size is large and there are many control genes. Interestingly, with few control genes the anti-conservative behavior gets worse as sample size increases.

7.5 Comparisons: Summary

Our empirical comparisons show that all methods for removing unwanted variation consistently improve on OLS in terms of ranking non-null genes vs null genes, with RUVB overall performing best here. In terms of calibration, several methods — including RUVB — were capable of providing good “typical” calibration across datasets. However, providing consistently correct calibration remains a challenge even in our relatively idealized scenarios. Use of control-gene calibration can be effective provided sample sizes are sufficiently large and sufficient control genes are available. In practice the main challenge here is likely to be specifying a sufficiently-large set of reliable control genes. Use of MAD calibration can ensure conservative behavior, but at a potential substantial loss in power to detect real effects.

8 Discussion

In this paper we developed a framework, RUV*, that both unites and generalizes different approaches to removing unwanted variation that use control genes. This unifying framework, which is based on viewing the problem as a matrix imputation problem, helps clarify connections between existing methods. In particular we provide conditions under which two popular methods, RUV2 and RUV4, are equivalent. The RUV* framework also preserves one of the advantages of existing RUV approaches — their modularity — which facilitates the development of novel methods based on existing matrix imputation algorithms. At the same time RUV* extends the RUV methods to allow for incorporating uncertainty in estimated unwanted variation. We provide one illustration of this via RUVB, a version of RUV* based on Bayesian factor analysis. In realistic simulations based on real data we found that RUVB is competitive with existing methods in terms of both power and calibration, although we also highlighted the challenges of providing consistently reliable calibration among data sets.

The methods developed in this paper are implemented in the R package vicar available at

https://github.com/dcgerard/vicar.

This package contains functions that easily allow analysts to include their own factor analysis code or (in the case of RUVB) their own prior specifications into the RUV pipelines. Code and instructions for reproducing the empirical evaluations in Section 7 are available at https://github.com/dcgerard/ruvb_sims.

9 Acknowledgments

This work was supported by NIH grant HG02585 and by a grant from the Gordon and Betty Moore Foundation (Grant GBMF #4559). Much of the original code for the simulated dataset generation in Section 7 was based on implementations by Mengyin Lu, to whom we are indebted. We would also like to express our sincere thanks to the authors of the cate [Wang and Zhao, 2015], ruv [Gagnon-Bartsch, 2015], and sva [Leek et al., 2016] R packages for making their work accessible for utilization or comparison. We also thank the Genotype-Tissue Expression Consortium [GTEx Consortium, 2015] for making their data available for analysis.

Appendix A Appendix

A.1 Approximate Posterior Inference

As discussed in Section 5.3, if we could calculate $p(\tilde{\boldsymbol{Y}}_{2\bar{\mathcal{C}}}|\mathcal{Y}_{m})$, where $\mathcal{Y}_{m}\mathrel{\mathop{\mathchar 58\relax}}=\{\boldsymbol{Y}_{2\mathcal{C}},\boldsymbol{Y}_{3\mathcal{C}},\boldsymbol{Y}_{3\bar{\mathcal{C}}}\}$, then inference on $[\boldsymbol{\beta}_{2}|\boldsymbol{Y}_{2\bar{\mathcal{C}}},\mathcal{Y}_{m}]$ would be straightforward — at least in principal if not in practice. That is, suppose $h(\tilde{\boldsymbol{Y}}_{2\bar{\mathcal{C}}})\mathrel{\mathop{\mathchar 58\relax}}=p(\tilde{\boldsymbol{Y}}_{2\bar{\mathcal{C}}}|\mathcal{Y}_{m})$, then one would simply use the likelihood $h(\boldsymbol{Y}_{2\bar{\mathcal{C}}}-\boldsymbol{R}_{22}\boldsymbol{\beta}_{2})$ and a user-provided prior $g(\cdot)$ over $\boldsymbol{\beta}_{2}$ to calculate a posterior and return posterior quantities.

However, to reap the benefits of modularity, we describe a procedure to fit the overall model (97) in two discrete steps: A factor analysis step using (103) and a regression step using (104). We begin with estimating the unwanted variation. Specifically, we suppose that one first assumes the model (103) where $\tilde{\boldsymbol{Y}}_{2\bar{\mathcal{C}}}$ is unobserved. The error $\boldsymbol{E}$ can follow any model a researcher desires. Though, of course, the rotation leading to (97) was derived by assuming Gaussian errors with independent rows (Section 2.1) and the appropriateness of different error models should be examined before use. We make the relatively weak assumption that model (103) is fit using any Bayesian procedure that yields a proper posterior and that the researcher can obtain samples from the posterior distribution $[\tilde{\boldsymbol{Y}}_{2\bar{\mathcal{C}}}|\mathcal{Y}_{m}]$. Call these posterior draws $\tilde{\boldsymbol{Y}}_{2\bar{\mathcal{C}}}^{(1)},\ldots,\tilde{\boldsymbol{Y}}_{2\bar{\mathcal{C}}}^{(t)}$.

After estimating the unwanted variation, we have a regression step where we estimate $\boldsymbol{\beta}_{2}$ using (104). Suppose we have any user-provided prior density over $\boldsymbol{\beta}_{2}$, say $g(\cdot)$. In order to reap the benefits of modularity, we need to derive a Bayesian procedure for approximating the posterior over $\boldsymbol{\beta}_{2}$ using just the samples $\tilde{\boldsymbol{Y}}_{2\bar{\mathcal{C}}}^{(1)},\ldots,\tilde{\boldsymbol{Y}}_{2\bar{\mathcal{C}}}^{(t)}$ from the previous step. To do so, we let $\hat{\boldsymbol{\beta}}_{2}^{(i)}\mathrel{\mathop{\mathchar 58\relax}}=\boldsymbol{R}_{22}^{-1}(\boldsymbol{Y}_{2\bar{\mathcal{C}}}-\tilde{\boldsymbol{Y}}_{2\bar{\mathcal{C}}}^{(i)})$ and note that $\hat{\boldsymbol{\beta}}_{2}^{(1)},\ldots,\hat{\boldsymbol{\beta}}_{2}^{(t)}$ are actually draws from $[\boldsymbol{\beta}_{2}|\boldsymbol{Y}_{2\bar{\mathcal{C}}},\mathcal{Y}_{m}]$ when using an (improper) uniform prior over $\boldsymbol{\beta}_{2}$. This follows because (104) is a location family. We can then weight these samples by the prior information $g(\hat{\boldsymbol{\beta}}_{2}^{(i)})$ to obtain draws from the posterior $[\boldsymbol{\beta}_{2}|\boldsymbol{Y}_{2\bar{\mathcal{C}}},\mathcal{Y}_{m}]$ when using $g(\cdot)$ as our prior density. This strategy of weighting samples from one distribution to approximate quantities from another distribution was discussed in Trotter and Tukey [1956]. What this means in practice is that given any function of $\boldsymbol{\beta}_{2}$, say $f(\cdot)$, we can approximate its posterior expectation consistently in the number of posterior draws, $t$, from the first step. That is,

[TABLE]

Some example calculations of useful posterior quantities are provided in Appendix A.5. We have given intuitive arguments here; formal arguments are given in Appendix A.2. A technical condition required for this approach to work is that $g(\cdot)$ be absolutely continuous with respect to the distribution of $[\boldsymbol{R}_{22}^{-1}(\boldsymbol{Y}_{2\bar{\mathcal{C}}}-\tilde{\boldsymbol{Y}}_{2\bar{\mathcal{C}}})|\mathcal{Y}_{m}]$. In the case when the errors $\boldsymbol{E}$ are Gaussian, it suffices to consider priors that are absolutely continuous with respect to Lebesgue measure.

Importantly, this two-step approach, though modular, is actually fitting the full model (97) as if done in one step. That is, nothing is lost in taking this two-step approach, except perhaps we could have found more efficient posterior approximations if the procedure was fit in one step. However, our approach is a contrast to RUV2, RUV3, and RUV4 which do not propagate the uncertainty in estimating the unwanted variation. This allows us to give more accurate quantities of uncertainty (Section 7).

To implement this approach in practice, we need to specify both a specific model for the unwanted variation (103) and a prior for $\boldsymbol{\beta}_{2}$. As a proof of concept we use a very simple Bayesian factor model with Gaussian errors (Appendix A.3) and an improper uniform prior on $\boldsymbol{\beta}_{2}$, which yields a proper proper posterior no matter the model for the unwanted variation (Appendix A.4). We note that although our model for the unwanted variation is based on a factor model, RUVB is neither a version of RUV4 nor RUV2.

A.2 Justification for Approximate Posterior Inference

In this section, we prove a general result that given a location family, we can approximate posterior expectations to any arbitrary level of precision using samples from the error distribution. We then connect this to the posterior approximation discussed in Appendix A.1. For data $\boldsymbol{y}\in\mathbb{R}^{d}$, suppose the model is

[TABLE]

where $h\mathrel{\mathop{\mathchar 58\relax}}\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ is bijective. Let $\boldsymbol{J}(\boldsymbol{z})$ be the Jacobian matrix of $h$. That is

[TABLE]

and let $|\boldsymbol{J}(\boldsymbol{z})|$ denote its determinant. Let $g$ be the prior of $\boldsymbol{\theta}$, which we assume is absolutely continuous with respect to the density of $h^{-1}(\boldsymbol{y}-\boldsymbol{e})$.

Theorem 3**.**

Let $\boldsymbol{e}_{1},\ldots,\boldsymbol{e}_{K}$ be i.i.d. random variables equal in distribution to $\boldsymbol{e}$. Let $u\mathrel{\mathop{\mathchar 58\relax}}\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ be a function. Let

[TABLE]

then

[TABLE]

Proof.

Let $f$ be the density of $\boldsymbol{e}$. Then $p(\boldsymbol{y}|\boldsymbol{\theta})=f(\boldsymbol{y}-h(\boldsymbol{\theta}))$ since $\boldsymbol{y}$ belongs to a location family. We have

[TABLE]

Line (126) follows by two applications of the weak law of large numbers followed by Slutsky’s theorem. Line (128) follows by a change of variables $\boldsymbol{e}=\boldsymbol{y}-h(\boldsymbol{z})$, the Jacobian of which is just $\boldsymbol{J}(\boldsymbol{z})$. The condition on the prior $g$ is used in (129) when we start considering $\boldsymbol{z}$ as a dummy variable for $\boldsymbol{\theta}$. To think intuitively about this condition on the prior, if the measure for $\boldsymbol{\theta}$ is non-zero on a set $\mathcal{A}$ in the parameterspace where the likelihood is non-zero but the measure is zero, then this approximation procedure would never sample $h^{-1}(\boldsymbol{y}-\boldsymbol{e}_{k})\in\mathcal{A}$. This is even though $\mathcal{A}$ does have a non-zero posterior probability. The absolute continuity condition prohibits this behavior. ∎

We now connect this general result to Appendix A.1. The $\boldsymbol{y}$, $\boldsymbol{\theta}$, and $\boldsymbol{e}$ in this section are the $\boldsymbol{Y}_{2\bar{\mathcal{C}}}$, $\boldsymbol{\beta}_{2}$, and $\tilde{\boldsymbol{Y}}_{2\bar{\mathcal{C}}}$, respectively, in Appendix A.1. So instead of having draws $\boldsymbol{e}_{1},\ldots,\boldsymbol{e}_{K}$, we have that $\tilde{\boldsymbol{Y}}_{2\bar{\mathcal{C}}}^{(1)},\ldots,\tilde{\boldsymbol{Y}}_{2\bar{\mathcal{C}}}^{(t)}$ are draws from $[\tilde{\boldsymbol{Y}}_{2\bar{\mathcal{C}}}|\mathcal{Y}_{m}]$. We also have that $h(\boldsymbol{\beta}_{2})=\boldsymbol{R}_{22}\boldsymbol{\beta}_{2}$, and so the determinant of the Jacobian is merely $|\boldsymbol{R}_{22}|^{p}$. Since this is a constant independent of $\boldsymbol{\beta}_{2}$, the Jacobians in the numerator and denominator cancel in (123).

A.3 Simple Bayesian Factor Analysis

In this section, we present a simple Bayesian factor analysis which we used in our implementation of RUVB. The factor model is

[TABLE]

We use conditionally conjugate priors on all parameters:

[TABLE]

where all hyper-parameters with a [math] subscript are assumed known. We let $\text{Gamma}(a,b)$ denote the Gamma distribution with mean $a/b$ and variance $a/b^{2}$. In the empirical evaluations of Section 7, we set the prior “sample sizes” to be small ($\rho_{0}=\alpha_{0}=0.1$) so that the prior is only weakly informative. The prior mean for the precisions ($\beta_{0}$) was set arbitrarily to 1. Following Ghosh and Dunson [2009], we set $\eta_{0}=\tau_{0}=1$.

The prior we use is similar in flavor to that of Ghosh and Dunson [2009], and like them we use the parameter expansion from [Gelman, 2006], which improves MCMC mixing. However, there are some important differences between our formulation and that of Ghosh and Dunson [2009]. First, we chose not to impose the usual identifiability conditions on the $\boldsymbol{F}$ matrix as we are not interested in the actual factors or factor loadings. Rather, our specifications induce a prior over the joint matrix of interest $\boldsymbol{L}\boldsymbol{F}$, which is identified. Second, the prior specification in Ghosh and Dunson [2009] is not order invariant. That is, their prior is influenced by the arbitrary ordering of the columns of $\boldsymbol{Y}$. This is a known problem [Leung and Drton, 2016] and we circumvent it by specifying an order invariant prior. Third, our prior specifications include hierarchical moderation of the variances, an important and well-established strategy in gene-expression studies [Smyth, 2004]. Finally, we chose to link the variances of the genes with those of the factors. This is approximately modeling a mean-variance relationship and we have found it to work well in practice. We emphasize here that we do not actually know the form of the unwanted variation in the simulations in Section 7, and so the good performance of RUVB there is not due to some “unfair advantage” in the choice of prior.

One step of a Gibbs sampler is presented in Procedure 5. Repeated applications of the steps in Procedure 5 will result in a Markov chain whose stationary distribution is the posterior distribution of the parameters in our model. As the calculations of the full conditionals for all parameters are fairly standard, we omit the detailed derivations.

A.4 Propriety of Posterior

In this section, we consider the propriety of the posterior $\boldsymbol{\beta}_{2}$ from Appendix A.1. To prove that the posterior of $\boldsymbol{\beta}_{2}$ under a uniform prior is proper, it suffices to consider the model (104), which we repeat here:

[TABLE]

Letting $\boldsymbol{A}\mathrel{\mathop{\mathchar 58\relax}}=\boldsymbol{R}_{22}^{-1}\boldsymbol{Y}_{\bar{\mathcal{C}}}$ and $\boldsymbol{C}=\boldsymbol{R}_{22}^{-1}\tilde{\boldsymbol{Y}}_{2\bar{\mathcal{C}}}$, we have

[TABLE]

Equation (140) represents a general location family with location parameter $\boldsymbol{\beta}_{2}$. It is not difficult to prove that using a uniform prior on the location parameter in a location family always results in a proper posterior.

Theorem 4**.**

Let $\boldsymbol{C}$ be a random variable with density $f$. Let $\boldsymbol{A}=\boldsymbol{\beta}_{2}+\boldsymbol{C}$. Suppose we place a uniform prior on $\boldsymbol{\beta}_{2}$. Then the posterior of $\boldsymbol{\beta}_{2}$ is proper.

Proof.

We first note that

[TABLE]

The density of $\boldsymbol{A}$ is just $f(\boldsymbol{A}-\boldsymbol{\beta}_{2})$. Since the prior of $\boldsymbol{\beta}_{2}$ is uniform, this means that the posterior of $\boldsymbol{\beta}_{2}$ given $\boldsymbol{A}$ is proportional to the likelihood $f(\boldsymbol{A}-\boldsymbol{\beta}_{2})$. Hence, making a change of variables of $\boldsymbol{\eta}=\boldsymbol{A}-\boldsymbol{\beta}_{2}$, we have

[TABLE]

∎

A.5 Posterior Summaries

Using (120), one can obtain approximations to posterior summaries quite easily. Let $\hat{\boldsymbol{\beta}}_{2}^{(i)}\mathrel{\mathop{\mathchar 58\relax}}=\boldsymbol{R}_{22}^{-1}(\boldsymbol{Y}_{2\bar{\mathcal{C}}}-\tilde{\boldsymbol{Y}}_{2\bar{\mathcal{C}}}^{(i)})$. The posterior mean may be approximated by

[TABLE]

Suppose one desires a $(1-\alpha)$ credible interval for $\beta_{2jk}$ for some $0<\alpha<1$. Sort the $\hat{\beta}_{2jk}^{(i)}$’s such that $\hat{\beta}_{2jk}^{(1)}<\hat{\beta}_{2jk}^{(2)}<\cdots<\hat{\beta}_{2jk}^{(t)}$. A $(1-\alpha)$ credible interval may be approximated by finding the $\ell,m\in\{1,\ldots,t\}$ such that

[TABLE]

then setting the $(1-\alpha)$ interval as $(\hat{\beta}_{2jk}^{(\ell)},\hat{\beta}_{2jk}^{(m)})$. Note that in (144) we have assumed that the $\hat{\boldsymbol{\beta}}_{2}^{(i)}$’s all have completely distinct elements. If the error distribution is Gaussian then we may do this without loss of generality as $\tilde{\boldsymbol{Y}}_{2\bar{\mathcal{C}}}^{(i)}$ is drawn from some convolution with a normal, and so is absolutely continuous with respect to Lebesgue measure.

Local false sign rates (lfsr’s) [Stephens, 2016] have recently been proposed as a way to measure the confidence in the sign of each effect. The intuition is that the lfsr is the probability of making an error if one makes their best guess about the sign of a parameter. Let

[TABLE]

Then the lfsr’s may be approximated by

[TABLE]

which simplifies under a uniform prior to

[TABLE]

Though, in many cases in Section 7, the Markov chain during the Bayesian factor analysis did not sample $\beta_{2jk}$’s of opposite sign. The estimate for the lfsr using (146) would then be 0. As it is often desirable to obtain a ranking of the most significant genes, this is an unappealing feature. We instead use a normal approximation to estimate the the lfsr’s, using the posterior means and standard deviations from the samples of the $\beta_{2jk}$’s.

A.6 Additional Simulation Considerations

To implement the methods in the simulation studies in Section 7, there are additional technicalities to be considered. Here, we briefly list out our choices.

The number of factors, $q$, is required to be known for all methods that we examined. There are many approaches to estimate this in the literature [for a review, see Owen and Wang, 2016]. We chose to use the parallel analysis approach of [Buja and Eyuboglu, 1992] as implemented in the num.sv function in the R package sva [Leek et al., 2016]. An alternative choice could have been the bi-cross-validation approach described in [Owen and Wang, 2016] and implemented in the cate R package [Wang and Zhao, 2015].

RUV2, RUV3, RUV4, and CATE all require specifying a factor analysis (Definition 1) for their respective first steps. For all methods, we used the truncated SVD (56)-(58). This was to make the methods more comparable. Though we note that Wang et al. [2015] also suggest using the quasi-maximum likelihood approach of Bai and Li [2012].

For RUVB, we ran each Gibbs sampler (Algorithm 5) for 12,500 iterations, dropping the first 2,500 iterations as burn-in. We kept every 10th sample to retain 1000 posterior draws from which we calculated the posterior summaries of Appendix A.5. Convergence diagnostics and other checks were implemented on a sample of the Markov chains in our simulation studies. We detected no problems with convergence (data not shown).

Since we explored many combinations of methods, we adopt the following notation:

• •

o = original variance estimates,

• •

m = MAD variance calibration,

• •

c = control-gene variance calibration,

• •

l = limma-moderated variances (EBVM),

• •

lb = limma-moderated variances (EBVM) before GLS (for either CATE or RUV3),

• •

la = limma-moderated variances (EBVM) after GLS (for either CATE or RUV3),

• •

d = delta-adjustment from CATE package (additive variance inflation),

• •

n = $t$ approximation for the likelihood in RUVB,

• •

nn = normal approximation for the likelihood in RUVB.

When comparing the AUC of different methods, certain combinations of methods theoretically have the same AUC. Specifically, applying MAD variance calibration (m) or control-gene variance calibration (c) does not alter the AUC of a method. Thus, we only need to compare one method of each of these groups to obtain comprehensive results on AUC performance. The members of these groups are those shown in Supplementary Figure 1.

The effect of library size (the total number of gene-reads in a sample) is a well-known source of bias in RNA-seq data [Dillies et al., 2013]. We do not explicitly adjust for library size. In this paragraph, we briefly argue that RUV-like methods can effectively account for library size. The usual pipeline to adjust for library size is to choose a constant for each sample, $c_{i}$, and divide each count in a sample by $c_{i}$. Many proposals have been made to estimate $c_{i}$ [Anders and Huber, 2010, Bullard et al., 2010, Robinson and Oshlack, 2010]. There are also variations on this pipeline. For example, others choose a constant for each sample, $c_{i}$, and include the $\boldsymbol{c}=(c_{1},\ldots,c_{n})^{\intercal}$ vector as a covariate in the regression model [Langmead et al., 2010]. In terms of our $\log_{2}$-count matrix of responses $\boldsymbol{Y}$, this corresponds to fitting the model

[TABLE]

where $\boldsymbol{c}$ is estimated independently of $\boldsymbol{X}$ in some (ad-hoc) fashion and $\boldsymbol{d}$ may or may not be assumed to be the vector of ones, $\boldsymbol{1}_{p}$. However, equation (148) is just a factor-augmented regression model, the same as (1). Methods that assume model (1) thus need not adjust for library size and need not choose one of the many procedures to estimate $\boldsymbol{c}$. That is, library size can be treated as just another source of unwanted variation.

A.7 Calibrated CATE

We can derive a multivariate version of (111) and formally justify it with maximum likelihood arguments via a modification of step 2 of Procedure 1. After step 1, we modify (10) and (11) to include a variance inflation parameter, $\lambda$.

[TABLE]

Step 2 of CATE simply calculates the MLE of $\boldsymbol{Z}_{2}$ in (10) and (11). Our calibration is simply finding the MLE’s of $\boldsymbol{Z}_{2}$ and $\lambda$ in (149) and (150), which can be done in closed form. The estimate of $\boldsymbol{Z}_{2}$ is unchanged (13). The MLE of $\lambda$ is

[TABLE]

Inference then proceeds by using $\hat{\lambda}\hat{\sigma}_{j}$ as the estimate of $\sigma_{j}$. Formula (151) is a multivariate version of equation (236) of Gagnon-Bartsch et al. [2013] (and (111)). Though while we derived (151) using an application of maximum likelihood, Gagnon-Bartsch et al. [2013] formulated their calibration using heuristic arguments.

Supplementary Materials