Power calculation for overall hypothesis testing with high-dimensional commensurate outcomes

1.1. Motivation

Sample size determination is one of the important aspects of study planning. An overestimation of sample size misuses research resources, while an underestimation results in lack of adequate statistical power. High-throughput methods used in metabolomic, genomic, and proteomic research generate high-dimensional data. When planning such studies, scientists increasingly emphasize pathway and set-based analysis, as opposed to simultaneous analysis of each individual metabolite, gene, and protein. Pathway analysis of functionally or structurally related molecules provides interpretable results that may not be easily accessible from simultaneous, individual testing. Goeman and Bühlmann [1] indicated that, in microarray analysis, the shift of the level of analysis from single genes to sets of related genes allows biologists to make use of previously accumulated biological knowledge for a more biology driven analysis.

Owing to the complexity of system biology, any metabolic, genetic, and proteomic pathways typically includes so many components (e.g., molecules) that statistical methods specialized for overall testing of high-dimensional outcomes are required. Concentrations or relative abundances of molecules are typically measured from one single high-throughput bioassay, leading to commensurate measurements sharing the same scale and units. An adequate account of the correlation between the high-dimensional commensurate outcomes can increase statistical power and in turn reduce the sample size required for a targeted power. Traditionally, multivariate methods were developed to account for such within-subject correlation (of any structure); however, the classical multivariate statistics become undefined when the number of variables exceeds sample size.

Many statistical methods have been proposed for a single overall test of high-dimensional differences in a biological pathway or gene set due to, for example, treatment and type of cell ([1–9], among others). However, most approaches apply to a narrow range of designs. Even fewer have power and sample size methods.

With high-dimensional outcomes, typical of genetic research, for example, sample size is often smaller than the number of variables because of cost considerations. Existing power and sample size methods for overall testing either are based on an assumption that sample size is greater than number of variables [10, 11] or rely on large sample Gaussian approximations [7, 8]. Neither approach may be suitable to allow samples with more variables than subjects. Moreover, the lack of free, user-friendly power and sample size software impedes study planning. In this paper, accurate power and sample size methods and software, applicable whether sample size is smaller or larger than the number of variables, will be developed. The new methods will be useful for planning studies for metabolic, genetic, and proteomic pathway analysis, with one-way or factorial between-subject design, adjustment for fixed covariates. Furthermore, the methods apply to any collection of fixed predictors in a general linear multivariate model, including incomplete designs. The methods allow computing power for any kind of general linear multivariate model test, including overall tests of mean vectors, within-subject contrasts, between-subject contrasts, and between-by-within interaction contrasts.

1.2. Previous related work

Driven by a collaboration to compare metabolic profiles before and after a vitamin B6 depleted diet, Chi et al. [5] proposed a new test that extends the existing ‘univariate approach’ to repeated measures (UNIREP) to analyzing high dimension, low sample size data (i.e., number of variables greater than sample size). On the basis of the framework of the general linear multivariate model, their test is flexible and can be applied to general designs, including parallel group comparison and evaluation of continuous covariates. Through a combination of analytic results and extensive simulations, they provided evidence that the new test has better control of type I error rate and power, when compared with methods proposed by Ahmad et al. [12], Srivastava and Du [7], and Srivastava and Fujikoshi [8]. Their method is implemented in a free, SAS/IML program called LINMOD, as well as being available in the SAS procedure PROC GLM (using uepsdef = CM option on the REPEATED statement).

Computer simulations provide a general tool for power analysis when the exact or an approximate null distribution of the statistic is available and the non-null distribution is unknown. However, simulations can introduce a substantial computational burden in high-dimensional settings. We derive the exact non-null distribution and a corresponding approximate distribution for the UNIREP statistic. The results provide fast and accurate power calculations. By the method of moments, our F approximation has properly adjusted degrees of freedom and reduces to the F approximation derived by Chi et al. [5] when the noncentrality parameter is zero. We describe the exact and approximate non-null distributions in terms of minimum sufficient sets of parameters, which specify a particular power computation. The theory assumes fixed predictors (i.e., with values known before data collection). As discussed by Glueck and Muller [13], different theory and computation of power analysis is needed when random predictors are present.

1.3. New sample size method and software

A minimum set of constants and parameters sufficient to specify the exact non-null distributions consists of (1) error degrees of freedom, (2) number of between-subject contrasts, (3) number of within-subject contrasts, (4) scaled variances of principal components of the hypothesis variables, and (5) multiple semi-partial correlations between principal components and the set of predictors tested in the hypothesis (adjusted for predictors in the model but not included in the test). We derive power-equivalent hypothesis testing scenarios that share the exact non-null distribution. We simulate only simple scenarios with canonical design matrices and diagonal covariance matrices. The analytic results on power equivalence guarantee that the conclusions from the simulations apply to scenarios with complex design matrices and any covariance structure.

The new power approximation uses only three parameters in lieu of the roughly p² needed for exact calculations. Simulation studies support the accuracy of the approximation. An updated version of SAS/IML software POWERLIB implements the new approximation in two ways. First, the program provides approximate power based on knowing the population parameters, including the entire covariance matrix and mean matrix, as is typical of power calculation software. Second, the program also allows using sample data (even with more variables than subjects) from a previous study to estimate the unknown parameters in the approximation and provide an estimated power for a new study.

We organize the rest of the paper as follows. In Section 2, we review the existing methods for testing general linear hypotheses and their distributional properties, especially with more variables than subjects. We then describe new exact and approximate results for the non-null distribution, highlight a number of power-equivalent hypothesis testing scenarios, and discuss power calculation when the population parameters are unknown. In Section 3, we provide a practical example for planning a study of metabolic consequences of vitamin B6 deficiency. In Section 4, we summarize results from extensive simulations for evaluating the accuracy of the approximations. In Section 5, we discuss implications of the new results and future research directions. We include all mathematical theorems and proofs in the Appendix.

2.1. Existing results for testing the general linear hypothesis

With rows of Y (N × p) corresponding to subjects (independent sampling units) and columns of Y corresponding to repeated measures or commensurate outcomes, the general linear multivariate model is given by

Here, X (N × q) is the design matrix, and B (q × p) is the primary parameters matrix. The Gaussian assumption leads to independently and identically distributed row_i (E)′ ~ 𝒩_p(0, Σ), equivalently E ~ 𝒩_N,p(0, I_N, Σ) for a matrix Gaussian distribution, as defined in [14, Chapter 8]. Table II summarizes the notation for the parameters and constants of the model and associated general linear hypothesis H₀ :Θ = Θ₀, with Θ = CBU the matrix of secondary parameters, C (a × q) the between-subject contrast matrix, and U (p × b) the within-subject contrast matrix.

A set of regularity conditions define estimable and testable Θ: (1) rank(C) = a ≤ q, (2) rank(U) = b ≤ p, and (3) C = C(X′ X)⁻(X′ X). Rows of the between-subject contrast matrix C define a contrasts in the predictor space, while columns of the within-subject contrast matrix U define b contrasts in the outcome space. In turn, the first two regularity conditions lead to defining an estimable and testable subspace. The third regularity condition is needed for the linear estimability of Θ ([15, p. 137] and [14, Section 11.4]). For an estimable Θ, both least squares and maximum likelihood estimation methods give Θ̂ = C B̃ U with B̃ = (X′ X)⁻X′ Y and (X′ X)⁻ the generalized inverse of X′ X. The general linear hypothesis H₀ : Θ = Θ₀ for testable Θ covers many designs for one, two, or many samples, including variations or analogs of repeated measures and multivariate analysis of variance, multivariate regression, discriminant analysis with two or more groups, and canonical correlation.

The UNIREP test statistic for the general linear multivariate model is proportional to the ratio of the trace of the hypothesis sum of squares matrix, Δ̂ = (Θ̂ − Θ₀) M⁻¹ (Θ̂ − Θ₀), to the trace of error sum of squares matrix, ν_eΣ̂_* = ν_eU′Σ̂U:

with Σ̂ = (Y − X B̃)′ (Y − X B̃) / ν_e and M = C (X′ X)⁻ C′. The statistic t_u is a one-to-one function of the sample estimate η̂_u = tr (Δ̂) = [tr (Δ̂) + ν_etr(Σ̂_*)] for the UNIREP measure of multivariate association:

Chi et al. [5] showed that under the null hypothesis H₀ : Θ = Θ₀, regardless of the ratio of sample size to the number of variables, the distribution function of t_u can be expressed exactly as

with y_kh ~ χ²(a) independent of y_ke ~ χ²(ν_e) for k ∈ {1, 2, ⋯, b}. With ${\pi }_k={\lambda }_k/\left({\sum }_{k=1}^b{\lambda }_k\right)$, the exact null distribution of t_u is invariant to global scale change in outcome space. From the prospective of principal component analysis, eigenvalues {λ_k} can be interpreted as variances of the principal components of hypothesis variables Y_u = YU, while the scaled eigenvalues {π_k} can be interpreted as proportions of variances explained by the principal components.

Chi et al. [5] extended the Box F approximation to the high dimension, low sample size setting:

Here, $\epsilon ={\left({\sum }_{k=1}^b{\lambda }_k\right)}^2/\left(b{\sum }_{k=1}^b{\lambda }_k^2\right)$ is the sphericity parameter quantifying the spread of population eigenvalues {λ_k} of Σ_*. They proposed a new estimator of ε, based on matching moments of the dual of the error covariance matrix. Their estimator is

with

Simulation results reported in Section 4 (including Tables III– VIII and Figure 2) of Chi et al. [5] demonstrate that their method (studied in the present paper) accurately controls type I error rate and has reasonable power even with a handful of subjects and thousands or more outcome variables. Chi et al. [5] reported empirical type I error rates for data simulated with b < ν_e (low dimension, high sample size) or b > ν_e (high dimension, low sample size). The range of designs included one-sample comparisons, two-sample comparisons, and three-sample comparisons with a covariate. Their test (T₂ in [5], also studied here) better controlled the type I error rate than the classic UNIREP tests when b > ν_e. The new test also outperformed the asymptotic tests proposed by Srivastava and Du [7] and Srivastava and Fujikoshi [8], especially with a small sample size.

2.2. New results: Exact non-null distribution

Analytic power calculation requires knowing the distribution of the statistic under the alternative hypothesis (e.g., H₁ : Θ ≠ Θ₀). In this section, we describe the exact noncentral distribution of t_u and its properties and delay the discussion about approximation to the next section.

Theorem 1 gives a simple expression for the exact noncentral distribution of t_u, which has a parallel form to the null distribution in Equation (5), except with y_kh ~ χ²(a, ω_k) for k ∈ {1, 2, ⋯, b}. The noncentral distribution is fully determined by specifying 3 constants and 2 · b parameters. The constants are determined by the design and hypothesis: (1) the error degrees of freedom, ν_e, (2) the number of between-subject contrasts, a, and (3) the number of within-subject contrasts, b. The unknown sets of sufficient parameters are {π_k} and {ω_k}.

The noncentrality parameters {ω_k} are the diagonal elements of Ω_t, the transformed, scaled non-centrality for the hypothesis sums of squares matrix Δ̂. Theorem 2(b) states that each noncentrality parameter is a one-to-one function of ${\rho }_k^2$, the squared multiple semi-partial correlation between principal components of hypothesis variables Y_u, and the set of predictors tested, with the predictors adjusted for all untested predictors in the model. The relationship is

Theorem 2(c) expresses the UNIREP measure of multivariate association in Equation (3) as

The squared multiple correlations $\{{\rho }_k^2\}$ are the univariate-multivariable coefficients of determination. They describe the proportion of variance in principal component k explained by the set of a predictors tested in the hypothesis (adjusted for predictors in the model but not included in the test). The interpretability of $\{{\rho }_k^2\}$ allows straightforward elicitation of {ω_k}. Together, {π_k} and $\{{\rho }_k^2\}$ define one set of sufficient parameters for exact noncentral distribution of t_u. They are the variance and correlation parameters characterizing the population structure in the principal component space. The sufficiency allows simplifying study planning and simulation designs because eliciting high-dimensional variance and covariance parameters can be a daunting task. The number of parameters greatly reduces from pq + p(p + 1)/2 in the original parameter space for B and Σ to 2 · b in the reduced space for {π_k} and $\{{\rho }_k^2\}$.

Any test that is based on a general linear multivariate model and its associated UNIREP statistic shares the same noncentral distribution and thus is power equivalent when N, ν_e, a, b, {π_k}, and {ω_k} are the same (Corollary 1). The equivalence holds if {λ_k} replaces {π_k} or $\{{\rho }_k^2\}$ replaces {ω_k}.

In Table II, three power-equivalent hypothesis testing scenarios are defined and depicted. Theorem 2(a) gives details of their statistical properties. The three scenarios share both {π_k} and {ω_k} in addition to constants N, ν_e, a, and b. Scenario S₁ has the original multivariate model defined in Equation (1), while scenarios S₂ and S₃ are based on transformed versions of the original model. Scenarios S2 has outcome variables Y₂ = YU ϒ = Y_u ϒ. With 𝓥(row_i (Y₂)) = Dg(λ), outcome variables for scenario S₂ are principal components of the hypothesis variables. The power equivalence between scenarios S₁ and S₂ allows assuming a diagonal hypothesis error covariance matrix Σ_* without loss of generality.

Following the notation in Table I, scenario S₃ in Table II has outcome variables Y₃ = T Y₂ with $\mathit{T}=\left[\begin{array}{cc}\hfill {\mathit{R}}_D^{\prime }\hfill & \hfill \mathbf{0}\hfill \\ \hfill \mathbf{0}\hfill & \hfill {\mathit{I}}_{{\nu }_e}\hfill \end{array}\right]{\mathit{L}}_X^{\prime }$ and R_D the eigenvectors of D′ D for D = CR_X+Dg(s_X+)⁻¹. Also, it has a simple, generic design matrix $\left[\begin{array}{c}\hfill {\mathit{I}}_r\hfill \\ \hfill \mathbf{0}\hfill \end{array}\right]$ and between-subject contrast matrix [I_a 0].With 𝓥(row_i(Y₃)) = Dg(λ), outcome variables for scenario S₃ are transformed components of the hypothesis variables. The power equivalence between scenarios S₁ and S₃ allows assuming a diagonal hypothesis error covariance matrix Σ_* and canonical design and contrast matrices without loss of generality. An important implication is that simulation results for diagonal Σ_* and canonical structure generalize to models with any covariance and design structure.

2.3. New results: Approximate non-null distribution

The lack of a closed form for the exact distribution of t_u leads to considering an approximation. By the method of moments, the distribution of t_u under the alternative hypothesis H₁ : Θ ≠ Θ₀ can be approximated by a noncentral F distribution, namely

Theorem 3 provides detailed expressions for parameters ε_n, ε_d = ε, and ω_u, all functions of {λ_k} and {ω_k}. Under the null hypothesis, ω_u = 0 and ε_n = ε, leading to the central F approximation in Equation (6). At the test size α, power can be approximated by computing 1 − Pr{F(abε_n, bν_eε_d, ω_u) ≤ f_crit} with Pr {F (abE (ε̃), ν_ebE(ε̃)) ≤ f_crit} = 1 − α. We extend the results in [11] to approximate E(ε̃) by $[({\nu }_a-2)({\nu }_a-4)/{\nu }_a^2]E({\mathrm{\epsilon ̃}}_{\mathrm{H}\mathrm{F}})$ with ε̃_HF = τ̂₁/ (bτ̂₂) = [ν_e + 1)bε̂ − 2] / [b (ν_e − b ε̂] and $\mathrm{\epsilon ̂}=b^{-1}{\text{tr}}^2\left({\mathbf{\Sigma ̂}}_{*}\right)/\text{tr}\left({\mathbf{\Sigma ̂}}_{*}^2\right)$.

The noncentral F approximation is fully determined by only three parameters, namely ε_n, ε_d = ε, and ω_u. Here, ω_u can be expressed as a function of ε_n and η_u, the UNIREP measure of multivariate association in Equation (12). Both ε_n and ω_u contains information pertaining to the magnitude and location of the association between the principal components of the hypothesis variables and the predictors tested. The minimum, sufficient parameter set for approximation (3 parameters) can be fully determined by the minimum, sufficient parameter set for exact power calculation (2b parameters), which in turn can be fully determined by the original parameter set (pq + p(p + 1)/2 parameters).

2.4. New results: Power calculation with estimated population parameters

As discussed in [16], statisticians frequently fix mean values and calculate power or sample size using a variance estimate from an existing study. For the general linear multivariate model given in Equation (1), fixed primary parameter matrix B leads to fixed secondary parameter matrix Θ and in turn to fixed hypothesis sum of squares matrix Δ. Adding a well-estimated Σ is sufficient for calculating an estimated power. A full rank estimate of Σ from an existing study allows estimating scaled eigenvalues of Σ_*, which in turn leads to specifications of the minimum sufficient parameters {π_k} and {ω_k} for the exact noncentral distribution. When the existing data have more variables than subjects, the estimator of Σ is singular and does not have adequate information to accurately provide estimates of the population eigenvalues [17].

In order to use data from an existing study with more variables than subjects, additional steps must be taken to estimate power for a new study. In practice, the three parameters {ε_d, ε_n, ω_u} suffice to specify the approximate noncentral F distribution (Equation (13)). By using Equations (8)–(10), we can first compute τ̂₁₀, τ̂₂₀, and ν_a0 from the existing data (hence the subscript 0) of sample size N₀, design matrix X₀, error degrees of freedom ν_e0 = N₀ − rank(X₀), and hypothesis error covariance matrix Σ̂_*0. We can then use Equations (7), (20), and (22) to estimate {ε_d, ε_n, ω_u} as

with Δ_t = (CBU − Θ₀)′ [C (X′_tX_t)⁻ C′]⁻¹ (CBU − Θ₀) and X_t the design matrix for the planned study. We note that the primary parameter matrix B is generally known and does not require estimation. The existing data are used only to estimate Σ by Σ̂₀, leading to which immediately gives estimates of Σ_* and functions of Σ_*, such as tr(Σ_*), tr²(Σ_*), and ε_d.With estimates of {ε_d, ε_n, ω_u} in hand, power can be approximated by computing 1 − Pr{F (abε̂_n, bν_eε̂_d, ω̂_u) ≤ f̂_crit} with Pr {F (abε̃, ν_ebε̃) ≤ f̂_crit} = 1 − α at the test size α. The form for the estimated critical value f̂_crit is based on the results from the extensive simulations in [18].

Variance estimates are known to have wide sampling distributions. Browne [19], Taylor and Muller [16, 20], and Muller and Pasour [21] described how to account for using estimates in power calculations in various settings with the univariate model. In the univariate setting, the 100(1 − β)th percentile of the variance estimates can be used to achieve the planned power in at least 100(1 − β)% of such studies. In the multivariate setting, which is the focus of the present manuscript, estimates for the b × b covariance matrix Σ_* lead to estimates of the parameters {ε_d, ε_n, ω_u} and the critical value f₀ (Equation (13)). Gribbin et al. [18] addressed the question for traditional low-dimensional designs with the sample size greater than the number of variables. Extending their results to the high dimension, low sample size design remains a topic for future research. Computing confidence bounds for power and sample size provides a natural solution to account for the variation and ensure achieving the planned power.

3.1. Design

Designing simulations with high-dimensional outcomes can be difficult as the number of variance and covariance parameters is of the magnitude of the squared number of outcomes, that is, b(b + 1)/2. On the basis of Corollary 2, power equivalence between scenarios S₁ and S₄ allows greatly simplifying the task. Scenario S₁ uses the general linear multivariate model Y = XB + E with no constraints on the mean of Y, namely XB, or on the covariance structure of the errors, namely Σ. In contrast, scenario S₄ is based on the model Y ₄ = X₄B₄ + E₄, which has four constraints: (1) independent errors, 𝓥[row_i (E′₄)] = Dg(π); (2) orthonormal predictors, X′₄X₄ = I_r ; (3) primary parameter matrix B₄ is sparse with all zeros except in the first row; and (4) simple contrast matrices, C₄ = [I_a 0] and U₄ = I_b. The power equivalence between scenarios S₁ and S₄ results from sharing a, b, ν_e, {π_k}, and {ω_k} (or $\{{\rho }_k^2\}$).

We designed our simulations on the basis of the data structures in scenario S₄. For both one-sample and two-sample comparisons, we adopted a five-way complete factorial design with factors b ∈ {64, 256, 1024}, N ∈ {10, 20, 40}, ε ∈ {0.27, 0.56, 0.76}, number of nonzero ${\rho }_k^2\in \{4,32\}$, and the location of nonzero ${\rho }_k^2$ at either the most dominant or middle components in hypothesis space. We also simulated a three-sample comparison with the same factorial design except for using N ∈ {12, 24, 48} with equal size groups. We varied the number of nonzero ${\rho }_k^2$ in order to address the power differences expected to result between concentrated and diffuse patterns of effect [22]. The location of nonzero ${\rho }_k^2$ dictates the sources of variation that are accounted for by predictors tested. The ordered, scaled variances π₁, …, π_b and τ were selected with π_k = (b − k + 1)^τ such that ε ∈ {0.27, 0.56, 0.76}, and ∑_k π_k = 1. The values of nonzero ${\rho }_k^2$ were set equal to each other and selected to achieve target power P ∈ {0.2, 0.5, 0.8, 0.9} for b = 64, and P ∈ {0.8, 0.9} for b ∈ {256, 1024}. All simulations used 100,000 replications and α = 0.05.

When b ∈ {256, 1024} and ε = 0.27, some of the variance ratios π_k/π₁ are numerical zeros, leading to an untestable hypothesis. The observation can be formalized as a generalization of Corollary 15.1 in [14] in the following way. A test for Θ₄ = C₄B₄U₄ is indistinguishable from a test for Θ_4* = C₄BU_4*, with U_4* = [I_b* 0]′ and b_* the number of the variance ratios π_k/π₁ greater than zeros, numerically. In order to assess whether the result was of practical concern, we compared power approximations for testing Θ₄ and Θ_4*. We observed that the approximated power values for the b-dimensional and b_*-dimensional tests coincide to at least the fourth decimal place. The expectation E(ε̃) varies at the third decimal place, leading to small differences in the critical values expected. We therefore report the simulation results for testing Θ₄ when the entire outcome space is considered.

3.2. One-sample problem

Table III displays a summary of absolute differences between the empirical and approximated power values. Performance of our proposed approximation was evaluated and summarized across a range of sample sizes N, covariance structures (governed by ε), and patterns of non-null effects (i.e., number and location of nonzero ${\rho }_k^2$). Across the board, our approximation produced absolute biases with both means and standard deviations less than 0.005. As the number of outcomes increased, the mean and standard deviation decreased. Overall, accuracy improved with number of outcomes, as well as sample size. The results are particularly useful for planning nutritional research when one-sample design is commonly used for studies comparing pre-intervention and post-intervention profiles.

3.3. Two-sample problem

Two-sample comparisons are very common in genomic and microarray research. Table IV summarizes the absolute bias when approximating power for testing differences between two groups. Both the means and standard deviations of the absolute biases were less than 0.01 when target power P ∈ {0.2, 0.5, 0.8, 0.9} and number of hypothesis outcomes b ∈ {64, 256, 1024}. Table V shows a detailed summary for b = 64 and P = 0.9 for a variety of sample sizes, noncentrality pattern, and covariance pattern. In general, as sample size increased, the bias decreased. The approximation worked slightly better when the group differences arose from the most dominant principal components of the hypothesis variables than from the middle dominant principal components. Overall, all documented biases were less than 0.03, giving approximated power in the range between 0.90 and 0.93.

3.4. Three-sample problem

Table VI summarizes the power accuracy for three-sample comparisons. Average bias was below 0.03, slightly higher than for one-sample and two-sample designs. Despite that, for power and sample size calculation, a bias at the level of 0.03 is usually acceptable.