 Methodology article
 Open Access
 Published:
Highly efficient hypothesis testing methods for regressiontype tests with correlated observations and heterogeneous variance structure
BMC Bioinformatics volume 20, Article number: 185 (2019)
Abstract
Background
For many practical hypothesis testing (HT) applications, the data are correlated and/or with heterogeneous variance structure. The regression ttest for weighted linear mixedeffects regression (LMER) is a legitimate choice because it accounts for complex covariance structure; however, high computational costs and occasional convergence issues make it impractical for analyzing highthroughput data. In this paper, we propose computationally efficient parametric and semiparametric tests based on a set of specialized matrix techniques dubbed as the PBtransformation. The PBtransformation has two advantages: 1. The PBtransformed data will have a scalar variancecovariance matrix. 2. The original HT problem will be reduced to an equivalent onesample HT problem. The transformed problem can then be approached by either the onesample Student’s ttest or Wilcoxon signed rank test.
Results
In simulation studies, the proposed methods outperform commonly used alternative methods under both normal and double exponential distributions. In particular, the PBtransformed ttest produces notably better results than the weighted LMER test, especially in the high correlation case, using only a small fraction of computational cost (3 versus 933 s). We apply these two methods to a set of RNAseq gene expression data collected in a breast cancer study. Pathway analyses show that the PBtransformed ttest reveals more biologically relevant findings in relation to breast cancer than the weighted LMER test.
Conclusions
As fast and numerically stable replacements for the weighted LMER test, the PBtransformed tests are especially suitable for “messy” highthroughput data that include both independent and matched/repeated samples. By using our method, the practitioners no longer have to choose between using partial data (applying paired tests to only the matched samples) or ignoring the correlation in the data (applying two sample tests to data with some correlated samples). Our method is implemented as an R package ‘PBtest’ and is available at https://github.com/yunzhang813/PBtestRPackage.
Background
Modern statistical applications are typically characterized by three major challenges: (a) highdimensionality; (b) heterogeneous variability of the data; and (c) correlation among observations. For example, numerous data sets are routinely produced by highthroughput technologies, such as microarray and nextgeneration sequencing, and it has become a common practice to investigate tens of thousands of hypotheses simultaneously for those data. When the classical i.i.d. assumption is met, the computational issue associated with highdimensional hypothesis testing (hereinafter, HT) problem is relatively easy to solve. As proof, R packages genefilter [1] and Rfast [2] implement vectorized computations of the Student’s and Welch’s ttests, respectively, both of which are hundreds times faster than the stock R function t.test(). However, it is common to observe heterogeneous variabilities between highthroughput samples, which violates the assumption of the Student’s ttest. For example, samples processed by a skillful technician usually have less variability than those processed by an inexperienced person. For twogroup comparisons, a special case of the heterogeneity of variance, i.e., samples in different groups have different variances, is well studied and commonly referred to as the BehrensFisher problem. The best known (approximate) parametric solution for this problem is the Welch’s ttest, which adjusts the degrees of freedom (hereinafter, DFs) associated with the tdistribution to compensate for the heteroscedasticity in the data. Unfortunately, the Welch’s ttest is not appropriate when the data have even more complicated variance structure. As an example, it is well known that the quality and variation of the RNAseq sample is largely affected by the total number of reads in the sequencing specimen [3, 4]. This quantity is also known as sequencing depth or library size, which may vary widely from sample to sample. Fortunately, such information is available a priori to data analyses. Several weighted methods [5–7] are proposed to utilize this information and make reliable statistical inference.
As the technology advances and the unit cost drops, immense amount of data are produced with even more complex variancecovariance structures. In multisite studies for big data consortium projects, investigators sometimes need to integrate omicsdata from different platforms (e.g. microarray or RNAseq for gene expression) and/or processed in different batches. Although many normalization [8–10] and batchcorrection methods [11–13] can be used to remove spurious bias, the heterogeneity of variance remains to be an issue. Besides, the clustering nature of these data may induce correlation among observations within one center/batch. Correlation may arise due to other reasons such as paired samples. For example, we downloaded a set of data for a comprehensive breast cancer study [14], which contain 226 samples including 153 tumor samples and 73 paired normal samples. Simple choices such as Welch’s ttest and paired ttest are not ideal for comparing the gene expression patterns between normal and cancerous samples, because they either ignore the correlations of the paired subjects or waste information contained in the unpaired subjects. To ignore the correlation and use a twosample test imprudently is harmful because it may increase the type I error rate extensively [15]. On the other hand, a paired test can only be applied to the matched samples, which almost certainly reduces the detection power. In general, data that involves two or more matched samples are called repeated measurements, and it is very common in practice to have some unmatched samples, also known as unbalanced study design.
One of the most versatile tools in statistics, the linear mixedeffects regression (LMER), provides an alternative inferential framework that accounts both unequal variances and certain practical correlation structures. The standard LMER can model the correlation by means of random effects. By adding weights to the model, the weighted LMER is able to capture very complex covariance structures in real applications. Although LMER has many nice theoretical properties, fitting it is computationally intensive. Currently, the best implementation is the R package lme4 [16], which is based on an iterative EM algorithm. For philosophical reasons, lme4 does not provide pvalues for the fitted models. The R package lmerTest [17] is the current practical standard to perform regression t and Ftests for lme4 outputs with appropriate DFs. A fast implementation of LMER is available in the Rfast package, which is based on highly optimized code in C++ [2]; however, this implementation does not allow for weights.
Many classical parametric tests, such as twosample and paired ttests, have their corresponding rankbased counterparts, i.e. the Wilcoxon ranksum test and the Wilcoxon signed rank test. A rankbased solution to the BehrensFisher problem can be derived based on the adaptive rank approach [18], but it was not designed for correlated observations. In recent years, researchers also extended rankbased tests to situations where both correlations and weights are presented. [19] derived the Wilcoxon ranksum statistic for correlated ranks, and [20] derived the weighted MannWithney U statistic for correlated data. These methods incorporate an interchangeable correlation in the whole dataset, and are less flexible for a combination of correlated and uncorrelated ranks. Lumley and Scott [21] proved the asymptotic properties for a class of weighted ranks under complex sampling, and pointed out that a reference tdistribution is more appropriate than the normal approximation for the Wilcoxon test when the design has low DFs. Their method is implemented in the svyranktest() function in R package survey. But most of the rankbased tests are designed for group comparisons; rankbased approaches for testing associations between two continuous variables with complex covariance structure are underdeveloped.
Based on a linear regression model, we propose two HT procedures (one parametric and one semiparametric) that utilize a priori information of the variance (weights) and correlation structure of the data. In “Methods” section, we design a linear map, dubbed as the “PBtransformation”, that a) transforms the original data with unequal variances and correlation into certain equivalent data that are independent and identically distributed; b) maps the original regressionlike HT problem into an equivalent onegroup testing problem. After the PBtransformation, classical parametric and rankbased tests with adjusted DFs are directly applicable. We also provide a moment estimator for the correlation coefficient for repeated measurements, which can be used to obtain an estimated covariance structure if it is not provided a priori. In “Simulations” section, we investigate the performance of the proposed methods using extensive simulations based on normal and double exponential distributions. We show that our methods have tighter control of type I error and more statistical power than a number of competing methods. In “A real data application” section, we apply the PBtransformed ttest to an RNAseq data for breast cancer. Utilizing the information of the paired samples and sequencing depths, our method selects more cancerspecific genes and fewer falsely significant genes (i.e. genes specific to other diseases) than the major competing method based on weighted LMER.
Lastly, computational efficiency is an important assessment of modern statistical methods. Depending on the number of hypotheses to be tested, our method can perform about 200 to 300 times faster than the weighted LMER approach in simulation studies and real data analyses. This efficiency makes our methods especially suitable for fast feature selection in highthroughput data analysis. We implement our methods in an R package called ’PBtest’, which is available at https://github.com/yunzhang813/PBtestRPackage.
Methods
Model framework
For clarity, we first present our main methodology development for a univariate regression problem. We will extend it to multiple regression problems in “Extension to multiple regressions” section.
Consider the following regressiontype HT problem:
Here, y is the response variable, x is the covariate, and ε is the error term that follows an ndimensional multivariate normal distribution \(\mathcal {N}\) with mean zero and a general variancecovariance matrix Σ. By considering a random variable Y in the ndimensional space, the above problem can also be stated as
In this model, μ is the intercept or grand mean that is a nuisance parameter, and β is the parameter of interest that quantifies the effect size. We express the variancecovariance matrix of ε in the form
where σ^{2} is a nonzero scalar that quantifies the magnitude of the covariance structure, and S is a symmetric, positivedefinite matrix that captures the shape of the covariance structure. Additional constraints are needed to determine σ^{2} and S; here, we choose a special form that can subsequently simplify our mathematical derivations. For any given Σ, define
From the above definition, we have the following nice property
Hereinafter, we refer to S the standardized structure matrix satisfying Eq. 5.
The proposed method
As a special case of Model (3), if S is proportional to I, the identity matrix, it is wellknown that regression ttest is a valid solution to this HT problem. If S≠I, e.g. the observed data are correlated and/or have heterogeneous variance structure, the assumptions of the standard ttest are violated. In this paper, we propose a linear transformation, namely \(\mathbf {P}\mathbf {B} : \mathbf {Y} \to \tilde {\mathbf {Y}}\), which transforms the original data to a new set of data that are independent and identically distributed. Furthermore, we prove that the transformed HT problem related to the new data is equivalent to the original problem, so that we can approach the original hypotheses using standard parametric (or later rankbased) tests with the new data.
To shed more lights on the proposed method, we first provide a graphical illustration in Fig. 1. The proposed procedure consists of three steps.

1
Estimate \(\hat \mu (\mathbf {Y})\) (i.e. the weighted mean of the original data), and subtract \(\hat \mu \) from all data. This process is an oblique (i.e. nonorthogonal) projection from \(\mathbb {R}^{n}\) to an (n−1)dimensional subspace of \(\mathbb {R}^{n}\). The intermediate data from this step is Y^{(1)} (i.e. the centered data). It’s clear that \(\mathbb {E} \mathbf {Y}^{(1)}\) is the origin of the reduced space if and only if H_{0} is true.

2
Use the eigendecomposition of the covariance matrix of Y^{(1)} to reshape its “elliptical” distribution to a “spherical” distribution. The intermediate data from this step is Y^{(2)}.

3
Use the QRdecomposition technique to find a unique rotation that transforms the original HT problem to an equivalent problem of testing for a constant deviation along the unit vector. The equivalent data generated from this step is \(\tilde {\mathbf {Y}}\), and the HT problem associated with \(\tilde {\mathbf {Y}}\) can be approached by existing parametric and rankbased methods.
In the proposed PBtransformation, Bmap performs both transformations in Step 1 and 2; Pmap from Step 3 is designed to improve the power of the proposed semiparametric test to be described in “A semiparametric generalization” section.
Centering data
Using weighted least squares, the mean estimation based on the original data is \(\hat {\mu } (\mathbf {Y})=\mathbf {1}' \mathbf {S}^{1} \mathbf {Y}\) (for details please see Additional file 1: Section S1.1). We subtract \(\hat {\mu }\) from all data points and define the centered data as
where J=1·1^{′} (i.e. a matrix of all 1’s). With some mathematical derivations (see Additional file 1: Section S1.1), we have
The Bmap
Now, we focus on S−J, which is the structure matrix of the centered data. Let TΛT^{′} denote the eigendecomposition of S−J. Since the data are centered, there are only n−1 nonzero eigenvalues. We express the decomposition as follows
where T_{n−1}∈M_{n×(n−1)} is a semiorthogonal matrix containing the first n−1 eigenvectors and Λ_{n−1}∈M_{(n−1)×(n−1)} is a diagonal matrix of nonzero eigenvalues. Based on Eq. 6, we define (see Additional file 1: Section S1.2)
so that \(\mathbf {Y}^{(2)} := \mathbf {B}\mathbf {Y} \in \mathbb {R}^{n1}\) have the following mean and covariance
We call the linear transformation represented by matrix B the “Bmap”. So far, we have centered the response variable, and standardized the general structure matrix S into the identity matrix I. However, the covariate and the alternative hypothesis in the original problem are also transformed by the Bmap. For normally distributed Y, the transformed HT problem in Eq. 7 is approachable by the regression ttest; however, there’s no appropriate rankbased counterpart. In order to conduct a rankbased test for Y with broader types of distribution, we propose the next transformation.
The Pmap
From Eq. 7, define the transformed covariate
We aim to find an orthogonal transformation that aligns z to 1_{n−1} in the reduced space. We construct such a transformation through the QR decomposition of the following object
where A∈M_{(n−1)×2} is a columnwise concatenation of vector z and the target vector 1_{n−1},Q∈M_{(n−1)×2} is a semiorthogonal matrix, and R∈M_{2×2} is an upper triangular matrix. We also define the following rotation matrix
Geometrically speaking, ξ= cosθ, where θ is the angle between z and 1_{n−1}.
With the above preparations, we have the following result.
Theorem 1
Matrix P:=I−QQ^{′}+Q Rot Q^{′}=I_{(n−1)×(n−1)}−Q(I_{2×2}−Rot)Q^{′} is the unique orthogonal transformation that satisfies the following properties:
Proof
See Additional file 1: Section 1.3. □
We call the linear transformation P defined by Theorem 1 the “Pmap”. Equation 9 ensures that this map is an orthogonal transformation. Equation 10 shows that the vector z is mapped to 1_{n−1} scaled by a factor ζ. Equation 11 is an invariant property in the linear subspace \(L_{\mathbf {z}}^{\perp }\), which is the orthogonal complement of the linear subspace spanned by 1_{n−1} and z, i.e. L_{z}=span(1_{n−1},z). This property defines a uniqueminimum map that only transforms the components of data in L_{z} and leaves the components in \(L_{\mathbf {z}}^{\perp }\) invariant. A similar idea of constructing rotation matrices has been used in [22].
With both B and P, we define the final transformed data as \( \tilde {\mathbf {Y}} := \mathbf {P} \mathbf {Y}^{(2)} = \mathbf {P} \mathbf {B} \mathbf {Y}\), which has the following joint distribution
The normality assumption implies that each \(\tilde Y_{i}\) follows an i.i.d. normal distribution, for i=1,⋯,n−1. The location parameter of the common marginal distribution is to be tested with unknown σ^{2}. Therefore, we can approach this equivalent HT problem with the classical onesample ttest and Wilcoxon signed rank test (more in “A semiparametric generalization” section).
Correlation estimation for repeated measurements
If Σ is unknown, we can decompose Σ in the following way
where W is a diagonal weight matrix and Cor is the corresponding correlation matrix. By definition, the weights are inversely proportional to the variance of the observations. In many real world applications including RNAseq analysis, those weights can be assigned a priori based on the quality of samples; but the correlation matrix Cor needs to be estimated from the data. In this section, we provide a momentbased estimator of Cor for a class of correlation structure that is commonly used for repeated measurements. This estimator does not require computationally intensive iterative algorithms.
Let Y be a collection of repeated measures from L subjects such that the observations from different subjects are independent. With an appropriate data rearrangement, the correlation matrix of Y can be written as a blockdiagonal matrix
We assume that the magnitude of correlation is the same across all blocks, and denote it by ρ. Each block can be expressed as \(\phantom {\dot {i}\!}\text {\texttt {Cor}}_{l}(\rho) = (1\rho) \mathbf {I}_{n_{l} \times n_{l}} + \rho \mathbf {J}_{n_{l} \times n_{l}}, \quad \text {for} \quad l = 1, \cdots, L,\) where n_{l} is the size of the lth block and \(n={\sum \nolimits }_{l=1}^{L} n_{l}\).
We estimate the correlation based on the weighted regression residuals \(\hat {\boldsymbol {\epsilon }}\) defined by Eq. (S3) in Additional file 1: Section S2.1. Define two forms of residual sum of squares
where \(\hat {\boldsymbol {\epsilon }}_{l}\) is the corresponding weighted residuals for the lth block. With these notations, we have the following Proposition.
Proposition 1
Denote \(\Sigma _{\epsilon } = \text {cov}(\hat {\boldsymbol {\epsilon }})\)and assume that for some nonzero σ^{2},
An estimator of ρ based on the first moments of SS_{1} and SS_{2} is
Moreover, if \(\hat {\boldsymbol {\epsilon }} \sim \mathcal {N}(\mathbf {0}, \Sigma _{\epsilon })\) and n_{1}=⋯=n_{L}=n/L (i.e. balanced design), the above estimator coincides with the maximum likelihood estimator of ρ, which has the form
Proof
See Additional file 1: Section S2.1. □
Standard correlation estimates are known to have downward bias [23], which can be corrected by the Olkin and Pratt’s method [24]. With this correction, our final correlation estimator is
Kenwardroger approximation to the degrees of freedom
The degree of freedom (DF) can have nontrivial impact on hypothesis testing when sample size is relatively small. Intuitively, a correlated observation carries “less information” than that of an independent observation. In such case, the effective DF is smaller than the apparent sample size. Simple examples include the twosample ttest and the paired ttest. Suppose there are n observations in each group, the former test has DF=2n−2 for i.i.d. observations, and the latter only has DF=n−1 because the observations are perfectly paired. These trivial examples indicate that we need to adjust the DF according to the correlation structure in our testing procedures.
We adopt the degrees of freedom approximation proposed by [25] (KR approximation henceforth) for the proposed tests. The KR approximation is a fast momentmatching method, which is efficiently implemented in R package pbkrtest[26]. In broad terms, we use the DF approximation as a tool to adjust the effective sample size when partially paired data are observed.
Alternative approach using mixedeffects model
As we mentioned in “Background” section, the HT problem stated in Model (3) for repeated measurements can also be approached by the linear mixedeffects regression (LMER) model. Suppose the ith observation is from the lth subject, we may fit the data with a random intercept model such that
where 1_{l} is the indicator function of the lth subject, \(\gamma \sim N\left (0, \sigma ^{2}_{\gamma }\right)\), and \(\epsilon _{i} \stackrel {i.i.d.}{\sim } N\left (0, \sigma ^{2}_{\epsilon }\right)\). The correlation is modeled as
The LMER model is typically fitted by a likelihood approach based on the EM algorithm. Weights can be incorporated in the likelihood function. The lmer() function in R package lme4 [16] provides a reference implementation for fitting the LMER model. The algorithm is an iterative procedure until convergence. Due to relatively high computational cost, the mixedeffects model has limited application in highthroughput data.
The R package lmerTest [17] performs hypothesis tests for lmer() outputs. By default, it adjusts the DF using the Satterthwaite’s approximation [27], and can optionally use the KR approximation.
A semiparametric generalization
In the above sections, we develop the PBtransformed ttest using linear algebra techniques. These techniques can be applied to nonnormal distributions to transform their mean vectors and covariance matrices as well. With the following proposition, we may extend the proposed method to an appropriate semiparametric distribution family. By considering the uncorrelated observations with equal variance as a second order approximation of the data that we are approaching, we can apply a rankbased test on the transformed data to test the original hypotheses. We call this procedure the PBtransformed Wilcoxon test.
Proposition 2
Let \(\check {\mathbf {Y}} := \large \left \{ \check {Y}_{1}, \dots, \check {Y}_{n1}\large \right \}\) be a collection of i.i.d. random variables with a common symmetric density function g(y), g(−y)=g(y). Assume that \(\mathbb {E} \check {Y}_{1} = 0\), \(\text {var}(\check {Y}_{1}) = \sigma ^{2}\). Let Y^{∗} be a random number that is independent of \(\check {\mathbf {Y}}\) and has zero mean and variance σ^{2}. For every symmetric semidefinite \(\mathbf {S} \in \mathrm {M}_{n \times n}, \mathbf {x} \in \mathbb {R}^{n}\) and \(\mu, \beta \in \mathbb {R}\), there exists a linear transformation \(\mathbf {D}: \mathbb {R}^{n1} \to \mathbb {R}^{n}\) and constants u,v, such that
is an ndimensional random vector with
Furthermore, if we apply the PBtransformation to Y, the result is a sequence of (n−1) equal variance and uncorrelated random variables with zero mean if and only if β=0.
Proof
See Additional file 1: Section S1.4. □
The essence of this Proposition is that, starting with an i.i.d. sequence of random variables with a symmetric common p.d.f., we can use linear transformations to generate a family of distributions that is expressive enough to include a nonnormal distribution with an arbitrary covariance matrix and a mean vector specified by the effect to be tested. This distribution family is semiparametric because: a) the “shape” of the density function, g(y), has infinite degrees of freedom; b) the “transformation” (D, u, and v) has only finite parameters.
As mentioned before, applying both the B and Pmaps enables us to use the Wilcoxon signed rank test for the hypotheses with this semiparametric distribution family. This approach has better power than the test with only the Bmap as shown in “Simulations” section. Once the PBtransformed data are obtained, we calculate the Wilcoxon signed rank statistic and follow the testing approach in [21], which is to approximate the asymptotic distribution of the test statistic by a tdistribution with an adjusted DF. Note that Wilcoxon signed rank test is only valid when the underlying distribution is symmetric; therefore, the symmetry assumption in Proposition 2 is necessary. In summary, this PBtransformed Wilcoxon test provides an approximate test (up to the second order moment) for data that follow a flexible semiparametric distributional model.
Extension to multiple regressions
In this section, we present an extension of the proposed methods for the following multiple regression
Here the error term ε is assumed to have zero mean but does not need to have scalar covariance matrix. For example, ε can be the summation of random effects and measurement errors in a typical LMER model with a form specified in Eq. 4.
To test the significance of β_{k},k=1,…,p, we need to specify two regression models, the null and alternative models. Here the alternative model is just the full Model (16), and the null model is a regression model for which the covariate matrix is X_{−k}, which is constructed by removing the kth covariate (X_{k}) from X
Compared with the original univariate problem, we see that the nuisance covariates in the multiple regression case are X_{−k}β_{−k} instead of 1μ in Eq. 1. Consequently, we need to replace the centering step by regressing out the linear effects of X_{−k}
The new Btransformation is defined as the eigendecomposition of cov(E)=σ^{2}(S−X_{−k}X−k′). The Ptransformation is derived the same as before, but with the new B matrix.
Simulations
We design two simulation scenarios for this study: SIM1 for completely paired group comparison, and SIM2 for regressiontype test with a continuous covariate. For both scenarios we consider three underlying distributions (normal, double exponential, and logistic) and four correlation levels (ρ=0.2, ρ=0.4,ρ=0.6, and ρ=0.8). We compare the parametric and rankbased PBtransformed test with oracle and estimated correlation to an incomplete survey of alternative methods. Each scenario was repeated 20 times and the results of ρ=0.2 and 0.8 for normal and double exponential distributions are summarized in Figs. 2 and 3, and Tables 1 and 2. See Additional file 1, Section S3 for more details about the simulation design, additional results of ρ=0.4 and 0.6, and results for logistic distribution.
Figures 2 and 3 are ROC curves for SIM1 and SIM2, respectively. In all simulations, the proposed PBtransformed tests outperform the competing methods.
The PBtransformed ttest has almost identical performance with oracle or estimated ρ. Using the estimated ρ slightly lowers the ROC curve of the PBtransformed Wilcoxon test compared with the oracle curve, but it still has a large advantage over other tests. Within the parametric framework, the weighted LMER has the best performance among the competing methods. It achieves similar performance as our proposed parametric test when the correlation coefficient is small; however, its performance deteriorates when the correlation is large. Judging from the ROC curves, among the competing methods, the svyranktest() is the best rankbased test for the group comparison problem, primarily because it is capable of incorporating the correlation information. However, it fails to control the typeI error, as shown in Table 1.
Tables 1 and 2 summarize the typeI error rate and power at the 5% significance level for SIM1 and SIM2, respectively. Overall, the PBtransformed tests achieve the highest power in all simulations. In most cases, the proposed tests tend to be conservative in the control of typeI error; and replacing the oracle ρ by the estimated \(\hat {\rho }\) does not have significant impact on the performance of PBtransformed tests. The only caveat is the rankbased test for the regressionlike problem. Currently, there’s no appropriate method designed for this type of problem. When the oracle correlation coefficient is provided to the PBtransformed Wilcoxon test, it has tight control of type I error. With uncertainty in the estimated correlation coefficient, our PBtransformed Wilcoxon test may suffer from slightly inflated type I errors; but it is still more conservative than its competitors. Of note, other solutions, such as the naive ttest and rankbased tests, may have little or no power for correlated data, though they may not have the lowest ROC curve.
Computational cost and degrees of freedom
We record the system time for testing 2000 simulated hypotheses using our method and lmer(), since they are the most appropriate methods for the simulated data with the best statistical performance. Our method takes less than 0.3 s with given Σ, and less than 0.9 s with the estimation step; lmer() takes 182 s. We use a MacBook Pro equipped with 2.3 GHz Intel Core i7 processor and 8GB RAM (R platform: x86_64darwin15.6.0). Of note, lmer() may fail to converge occasionally, e.g. 0 – 25 failures (out of 2,000) in each repetition of our simulations. We resort to a try/catch structure in the R script to prevent these convergence issues from terminating the main loop.
We also check the degrees of freedom in all applicable tests. In this section, we report the DFs used/adjusted in SIM1, i.e. the completely paired group comparison. Recall that n=40 with n_{A}=n_{B}=20. It is straightforward to calculate the DFs used in the twosample ttest and the paired ttest, which are 38 and 19, respectively. Using lmerTest() (weighted LMER) with default parameters, it returns the mean DF = 35.51 with a large range (min = 4.77, max = 38) from the simulated data with ρ=0.2. Using the oracle Σ_{SIM}, our method returns the adjusted DF = 14.35; if the covariance matrix is estimated, our method returns the mean DF = 14.38 with high consistency (min = 14.36, max = 14.42). When ρ=0.8, the adjusted DFs become smaller. The weighted LMER returns the mean DF = 20.63 (min = 4.03, max = 38). Our method returns DF = 12.48 for the oracle covariance, and mean DF = 12.56 (min = 12.55, max = 12.57) for the estimated covariance. Also, the rankbased test svyranktest() returns a DF for its tdistribution approximation, which is 18 for both small and large correlations.
A real data application
We download a set of RNAseq gene expression data from The Cancer Genome Atlas (TCGA) [14] (see Additional file 1: Section S4). The data are sequenced on the Illumina GA platform with tissues collected from breast cancer subjects. In particular, we select 28 samples from the tissue source site “BH”, which are controlled for white female subjects with the HER2positive (HER2+) [28] biomarkers. After data preprocessing based on nonspecific filtering (see Additional file 1: Section S4.1), a total number of 11,453 genes are kept for subsequent analyses. Among these data are 10 pairs of matched tumor and normal samples, 6 unmatched tumor samples, and 2 unmatched normal samples. Using Eq. 13, the estimated correlation between matched samples across all genes is \(\hat {\rho } = 0.10\).
The sequencing depths of the selected samples range from 23.80 million reads to 76.08 million reads. As mentioned before, the more reads are sequenced, the better is the quality of RNAseq data [4]; thus it is reasonable to weigh samples by their sequencing depths. Since this quantity is typically measured in million reads, we set the weights
for i=1,⋯,28.
With the above correlation estimate and weights, we obtained the covariance structure using Eq. 12. For properly preprocessed sequencing data, a proximity of normality can be warranted [29]. We applied the PBtransformed ttest and the weighted LMER on the data.
Based on the simulations, we expect that if correlation is small, the PBtransformed ttest should have tighter control of false positives than alternative methods. At 5% false discovery rate (FDR) level combined with a foldchange (FC) criterion (FC<0.5 or FC>2), the PBtransformed ttest selected 3,340 DEGs and the weighted LMER selected 3,485 DEGs (for biological insights of the DEG lists, see Additional file 1: Section S4.4).
To make the comparison between these two methods more fair and meaningful, we focus on studying the biological annotations of the top 2,000 genes from each DEG list. Specifically, we apply the gene set analysis tool DAVID [30] to the 147 genes that uniquely belong to one list. Both Gene Ontology (GO) biological processes [31] and KEGG pathways [32] are used for functional annotations. Terms identified based on the 147 unique genes in each DEG list are recorded in Additional file 1: Table S6. We further pin down two gene lists, which consist of genes that participate in more than five annotation terms in the above table: there are 11 such genes (PIK3R2, AKT3, MAPK13, PDGFRA, ADCY3, SHC2, CXCL12, CXCR4, GAB2, GAS6, and MYL9) for the PBtransformed ttest, and six (COX6B1, HSPA5, COX4I2, COX5A, UQCR10, and ERN1) for the weighted LMER. Expression level of these genes are plotted in Fig. 4. These DEGs are biologically important because they are involved in multiple biological pathways/ontology terms.
Those 11 genes uniquely identified by the PBtransformed ttest are known to be involved in cell survival, proliferation and migration. The CXCR4CXCL12 chemokine signaling pathway is one of the deregulated signaling pathway uniquely identified by PBtransformed ttest in HER2+ breast cancer cells. This pathway is known to play a crucial role in promoting breast cancer metastasis and has been reported to be associated with poor prognosis [33, 34]. Compared with the stateoftheart method (weighted LMER), the PBtransformed ttest identifies more genes whose protein products can be targeted by pharmaceutical inhibitors. CXCR4 inhibitors have already demonstrated promising antitumor activities against breast [35, 36], prostrate [37] and lung [38] cancers. Additional downstream signaling molecules identified by our analysis to be significantly associated with HER2+ breast tumor such as PI3K, p38, adaptor molecule GAB2 and SHC2 can also be potential therapeutic targets for selectively eliminating cancer cells. Please refer to Additional file 1: Section S4.5 for full list of functional annotation terms.
Discussion
In this paper, we present a data transformation technique that can be used in conjunction with both the Student’s ttype test and rankbased test. In the simulation studies, our proposed tests outperform the classical tests (e.g. twosample/regreesion ttest and Wilcoxon ranksum test) by a large margin. In a sense, this superiority is expected, because the classical methods do not consider the correlation nor heteroscedasticity of the data.
In our opinion, the most practical comparison in this study is the one between the PBtransformed ttest and the weighted LMER. The fact that the PBtransformed ttest outperforms the weighted LMER, and this advantage is more pronounced for data with higher correlation (see e.g., Figs. 2 and 3), is the highlight of this study, which may have profound implications for applied statistical practice.
We believe the following reasons may explain the advantages of the PBtransformed tests. 1. As reported in “Computational cost and degrees of freedom” section, the default degrees of freedom approximation in lmerTest varies dramatically, as oppose to very stable degrees of freedom approximation in our method. 2. Our momentbased correlation estimator is better than the LMER correlation estimator (see Additional file 1: Section S2.2). One possible explanation is that LMER depends on nonlinear optimizer, which may not always converge to the global maximum likelihood. 3. In a minor way but related to 2, lmer() fails to converge to even a local maximum in certain rare cases.
Another major contribution of our method is that the transformationbased approach is computationally much more efficient than the EM algorithm used in LMER, which is an important advantage in highthroughput data analysis. Recall that in simulation studies, PBtransformed ttest is approximately 200 times faster than the weighted LMER approach. As an additional evidence, to test the 11,453 genes in the real data study, it takes 933 s using the weighted LMER, and only 3 s using our method, which is more than 300 times faster.
Nonetheless, we want to emphasize that, by no means, our method is a replacement for LMER. The mixedeffects model is a comprehensive statistical inference framework that includes parameter estimation, model fitting (and possibly model selection), hypothesis testing, among other things; whereas our methods are only designed for the hypothesis testing. We envision that in a typical highthroughput data application, an investigator may quickly run PBtransformed ttest to identify important features first, then apply lme4 to fit mixed effects models for those selected features. In this way, he/she enjoys both the computational efficiency of our method and the comprehensive results provided by a full LMER model.
In “Extension to multiple regressions” section, we extend the PBtransformed tests for multiple regressions. We must point out two weaknesses in this approach. 1. The proposed extension is comparable to the regression ttest for individual covariates, not the ANOVA Ftest for the significance of several covariates simultaneously. In fact, the Bmap can be defined in this case so we can define a transformed parametric test easily; but there is no clear counterpart for the Pmap, which is needed to overcome the identifiability issue for the semiparametric generalization. 2. The performance of PBtransformations depends on a good estimation of S, the shape of the covariance matrix of the observations. Currently, our momentbased estimator only works for problems with just one random intercept, which is only appropriate for relatively simple longitudinal experiments. It is a challenging problem to estimate the complex covariance structure for general LMER models (e.g., one random intercept plus several random slopes), and we think it can be a nice and ambitious research project for us in the near future.
Numerically, the PBtransformed ttest provides the same test statistic and degrees of freedom as those from the paired ttest for perfectly paired data and the regression ttest for i.i.d. data. In this sense, the PBtransformed ttest is a legitimate generalization of these two classical tests. The rankbased test is slightly different from the classical ones, since we used a tdistribution approximation instead of a normal approximation for the rankbased statistic. The tdistribution approximation is preferred for correlated data because the effective sample size may be small even in a large dataset [21].
Recall that the PBtransformation is designed in a way that the transformed data have the desired first and second order moments. For nonnormal distributions, the transformed samples may not have the same higher order moments. Note that, the Pmap is currently defined in part by Eq. (11), the minimum action principle. Without this constraint, we will have some extra freedom in choosing the Pmap. In the future development, we will consider using this extra freedom of orthogonal transformation to minimize the discrepancy of higher order moments of the transformed samples for the semiparametric distribution family. This would require an optimization procedure on a submanifold of the orthogonal group, which may be computationally expensive. The advantage is that, by making the higher order moments more homogeneous across the transformed data, we may be able to further improve the statistical performance of the PBtransformed Wilcoxon test.
In this study, we presented an example in RNAseq data analysis. In recent bioinformatics research, advanced methods such as normalization and batcheffect correction were developed to deal with data heterogeneities in bioassays. While most of these approaches are focused on the first moment (i.e. correction for bias in the mean values), our approach provides a different perspective based on the second order moments (i.e. the covariance structure). The dramatic computational efficiency boost of our method also opens the door for investigators to use the PBtransformed tests for ultrahighdimensional data analysis, such as longitudinal studies of diffusion tensor imaging data at the voxellevel [39–41], in which about one million hypotheses need to be tested simultaneously. Finally, we think the PBtransformed Wilcoxon test can also be used in metaanalysis to combine results from several studies with high betweensite variability and certain correlation structure due to, e.g., site and subjectspecific random effects.
Abbreviations
 HT:

Hypothesis testing
 LMER:

Linear mixed effects regression
 DF:

Degrees of freedom
 KR:

KenwardRoger approximation
 TCGA:

The Cancer Genome Atlas
 DAVID:

The Database for Annotation, Visualization and Integrated Discovery
 GO:

Gene ontology
 KEGG:

Kyoto encyclopedia of genes and genomes
 DEG:

Differential expressed genes
References
 1
Gentleman R, Carey V, Huber W, Hahne F. Genefilter: genefilter: methods for filtering genes from highthroughput experiments. R package version 1.60.0. 2017.
 2
Papadakis M, Tsagris M, Dimitriadis M, Tsamardinos I, Fasiolo M, Borboudakis G, Burkardt J. Rfast: Fast r functions. R package version, 1.5. 2017;1(5).
 3
Wang L, Wang S, Li W. Rseqc: quality control of rnaseq experiments. Bioinformatics. 2012; 28(16):2184–5.
 4
Sims D, Sudbery I, Ilott NE, Heger A, Ponting CP. Sequencing depth and coverage: key considerations in genomic analyses. Nat Rev Genet. 2014; 15(2):121–32.
 5
Law CW, Chen Y, Shi W, Smyth GK. Voom: precision weights unlock linear model analysis tools for rnaseq read counts. Genome Biol. 2014; 15(2):29.
 6
Zhou X, Lindsay H, Robinson MD. Robustly detecting differential expression in rna sequencing data using observation weights. Nucleic Acids Res. 2014; 42(11):91.
 7
Liu R, Holik AZ, Su S, Jansz N, Chen K, Leong HS, Blewitt ME, AsselinLabat ML, Smyth GK, Ritchie ME. Why weight? modelling sample and observational level variability improves power in rnaseq analyses. Nucleic Acids Res. 2015; 43(15):97.
 8
Robinson MD, Oshlack A. A scaling normalization method for differential expression analysis of rnaseq data. Genome Biol. 2010; 11(3):25.
 9
Risso D, Ngai J, Speed TP, Dudoit S. Normalization of rnaseq data using factor analysis of control genes or samples. Nat Biotechnol. 2014; 32(9):896–902.
 10
Liu Y, Zhang J, Qiu X. Superdelta: a new differential gene expression analysis procedure with robust data normalization. BMC Bioinformatics. 2017; 18(1):582.
 11
Johnson WE, Li C, Rabinovic A. Adjusting batch effects in microarray expression data using empirical bayes methods. Biostatistics. 2007; 8(1):118–27.
 12
Leek JT, Storey JD. Capturing heterogeneity in gene expression studies by surrogate variable analysis. PLoS Genet. 2007; 3(9):161.
 13
Hardcastle TJ, Kelly KA. bayseq: empirical bayesian methods for identifying differential expression in sequence count data. BMC Bioinformatics. 2010; 11(1):422.
 14
Cancer Genome Atlas Network T. Comprehensive molecular portraits of human breast tumors. Nature. 2012; 490(7418):61.
 15
Walsh JE. Concerning the effect of intraclass correlation on certain significance tests. Ann Math Stat. 1947; 18(1):88–96.
 16
Bates D, Mächler M, Bolker B, Walker S. Fitting linear mixedeffects models using lme4. arXiv preprint arXiv:1406.5823. 2014.
 17
Kuznetsova A, Brockhoff PB, Christensen RHB. lmerTest package: Tests in linear mixed effects models. J Stat Softw. 2017; 82(13):1–26. https://doi.org/10.18637/jss.v082.i13.
 18
Sidak Z, Sen PK, Hajek J. Theory of Rank Tests. San Diego: Academic press; 1999.
 19
Barry WT, Nobel AB, Wright FA. A statistical framework for testing functional categories in microarray data. Ann Appl Stat. 2008; 2(1):286–315.
 20
Zhang Y, Topham DJ, Thakar J, Qiu X. FunnelGSEA: Functional elasticnet regression in timecourse gene set enrichment analysis. Bioinformatics. 2017; 33(13):1944–52.
 21
Lumley T, Scott AJ. Twosample rank tests under complex sampling. Biometrika. 2013; 100(4):831–42.
 22
Amaral GA, Dryden I, Wood ATA. Pivotal bootstrap methods for ksample problems in directional statistics and shape analysis. J Am Stat Assoc. 2007; 102(478):695–707.
 23
Zimmerman DW, Zumbo BD, Williams RH. Bias in estimation and hypothesis testing of correlation. Psicológica. 2003; 24(1):133–159.
 24
Olkin I, Pratt JW. Unbiased estimation of certain correlation coefficients. Ann Math Stat. 1958; 29(1):201–211.
 25
Kenward MG, Roger JH. Small sample inference for fixed effects from restricted maximum likelihood. Biometrics. 1997; 53(3):983–997.
 26
Halekoh U, Hojsgaard S. A kenwardroger approximation and parametric bootstrap methods for tests in linear mixed models – the r package pbkrtest. J Stat Softw. 2014; 59(9):1–30.
 27
Satterthwaite FE. Synthesis of variance. Psychometrika. 1941; 6(5):309–16.
 28
Burstein HJ. The distinctive nature of her2positive breast cancers. N Engl J Med. 2005; 353(16):1652–4.
 29
Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK. limma powers differential expression analyses for rnasequencing and microarray studies. Nucleic Acids Res. 2015; 43(7):47.
 30
Huang DW, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using david bioinformatics resources. Nat Protoc. 2009; 4(1):44–57.
 31
Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, Davis AP, Dolinski K, Dwight SS, Eppig JT, et al.Gene ontology: tool for the unification of biology. Nat Genet. 2000; 25(1):25–9.
 32
Kanehisa M, Goto S. Kegg: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000; 28(1):27–30.
 33
Sun Y, Mao X, Fan C, Liu C, Guo A, Guan S, Jin Q, Li B, Yao F, Jin F. Cxcl12cxcr4 axis promotes the natural selection of breast cancer cell metastasis. Tumor Biol. 2014; 35(8):7765–73.
 34
Müller A, Homey B, Soto H, Ge N, Catron D, Buchanan ME, McClanahan T, Murphy E, Yuan W, Wagner SN, et al.Involvement of chemokine receptors in breast cancer metastasis. Nature. 2001; 410(6824):50.
 35
Huang EH, Singh B, Cristofanilli M, Gelovani J, Wei C, Vincent L, Cook KR, Lucci A. A cxcr4 antagonist ctce9908 inhibits primary tumor growth and metastasis of breast cancer1. J Surg Res. 2009; 155(2):231–6.
 36
Chittasupho C, Anuchapreeda S, Sarisuta N. Cxcr4 targeted dendrimer for anticancer drug delivery and breast cancer cell migration inhibition. Eur J Pharm Biopharm. 2017; 119:310–21.
 37
Wong D, Kandagatla P, Korz W, Chinni SR. Targeting cxcr4 with ctce9908 inhibits prostate tumor metastasis. BMC Urol. 2014; 14(1):12.
 38
Taromi S, Kayser G, Catusse J, von Elverfeldt D, Reichardt W, Braun F, Weber WA, Zeiser R, Burger M. Cxcr4 antagonists suppress small cell lung cancer progression. Oncotarget. 2016; 7(51):85185.
 39
Zhu T, Hu R, Tian W, Ekholm S, Schifitto G, Qiu X, Zhong J. Spatial regression analysis of diffusion tensor imaging (spread) for longitudinal progression of neurodegenerative disease in individual subjects. Magn Reson Imaging. 2013; 31(10):1657–67.
 40
Liu B, Qiu X, Zhu T, Tian W, Hu R, Ekholm S, Schifitto G, Zhong J. Improved spatial regression analysis of diffusion tensor imaging for lesion detection during longitudinal progression of multiple sclerosis in individual subjects. Phys Med Biol. 2016; 61(6):2497.
 41
Liu B, Qiu X, Zhu T, Tian W, Hu R, Ekholm S, Schifitto G, Zhong J. Spatial regression analysis of serial dti for subjectspecific longitudinal changes of neurodegenerative disease. NeuroImage Clin. 2016; 11:291–301.
Acknowledgments
Not applicable.
Funding
Research reported in this publication was supported in part by the National Institute of Environmental Health Sciences of the National Institutes of Health (NIH) under award number T32ES007271, the University of Rochester CTSA award number UL1 TR002001 from the National Center for Advancing Translational Sciences of the National Institutes of Health, the University of Rochester Center for AIDS Research (NIH 5 P30 AI07849808), and Respiratory Pathogens Research Center (NIAID contract number HHSN272201200005C). The content is solely the responsibility of the authors and does not necessarily represent the official views of the NIH.
Availability of data and materials
The methods are implemented in R package PBtest, freely and publicly available at https://github.com/yunzhang813/PBtestRPackage.
Author information
Affiliations
Contributions
XQ was responsible for the study design. YZ implemented the proposed method and performed data analysis. GB, DT, and AF provided biological interpretations of the real data. YZ. and XQ wrote the manuscript. All five authors revised the manuscript and approved the final version.
Corresponding author
Correspondence to Xing Qiu.
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Additional file
Additional file 1
This file contains: a) proofs of the main theorems; b) details of the momentbased correlation estimator; c) details of simulation design; and d) additional information about the real data analysis. (PDF 3523 kb)
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License(http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver(http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Cite this article
Zhang, Y., Bandyopadhyay, G., Topham, D.J. et al. Highly efficient hypothesis testing methods for regressiontype tests with correlated observations and heterogeneous variance structure. BMC Bioinformatics 20, 185 (2019). https://doi.org/10.1186/s1285901927838
Received:
Accepted:
Published:
Keywords
 Hypothesis testing
 Matrix decomposition
 Orthogonal transformation
 RNAseq
 Rotated test