Confidence regions for repeated measures ANOVA power curves based on estimated covariance.

BACKGROUND: Using covariance or mean estimates from previous data introduces randomness into each power value in a power curve. Creating confidence intervals about the power estimates improves study planning by allowing scientists to account for the uncertainty in the power estimates. Driving examples arise in many imaging applications.
METHODS: We use both analytical and Monte Carlo simulation methods. Our analytical derivations apply to power for tests with the univariate approach to repeated measures (UNIREP). Approximate confidence intervals and regions for power based on an estimated covariance matrix and fixed means are described. Extensive simulations are used to examine the properties of the approximations.
RESULTS: Closed-form expressions are given for approximate power and confidence intervals and regions. Monte Carlo simulations support the accuracy of the approximations for practical ranges of sample size, rank of the design matrix, error degrees of freedom, and the amount of deviation from sphericity. The new methods provide accurate coverage probabilities for all four UNIREP tests, even for small sample sizes. Accuracy is higher for higher power values than for lower power values, making the methods especially useful in practical research conditions. The new techniques allow the plotting of power confidence regions around an estimated power curve, an approach that has been well received by researchers. Free software makes the new methods readily available.
CONCLUSIONS: The new techniques allow a convenient way to account for the uncertainty of using an estimated covariance matrix in choosing a sample size for a repeated measures ANOVA design. Medical imaging and many other types of healthcare research often use repeated measures ANOVA.

Background

Motivation

Computing power for a linear model involving repeated measures requires specifying a set of means and a covariance matrix. Scientists usually feel comfortable specifying a pattern of means that corresponds to a difference of clinical or scientific importance. However, specifying plausible variance and covariance values usually requires estimates from a previous study.

Using data from a previous study to estimate the covariance matrix makes the power value a random variable. Kraemer, et al. [1] noted that if the estimated variance is too small, i.e., when the pilot study is overly favorable, power will be over-estimated. Conversely, if the estimated variance is too large (pilot study overly conservative), power will be under-estimated. Maxwell [2] conducted simulations to illustrate the amount of bias that can occur. Taylor and Muller [3] and Muller and Pasour [4] derived exact distributions of noncentrality and power in univariate linear models based on all combinations of estimated variance and means. The results account for power computed conditional on a previous result being significant, or conditional on a previous result being non-significant. The former creates optimistic bias, while the latter creates pessimistic bias.

Providing confidence intervals to account for the uncertainty inherent in the random power values would be useful for study planning. For example, a lower bound for power would allow stating that a test has power of at least “P” to detect an effect, with a specified confidence. A confidence region for a power curve would be even more informative.

Medical imaging research motivated the work here because it often generates the type of complete data that can be handled with the univariate approach to repeated measures (UNIREP). Muller, et al. [5] reviewed the advantages gained by being able to use the UNIREP model, a special case of the general linear mixed model. The same authors described accurate and convenient power approximations for UNIREP analysis. The four UNIREP tests, Box conservative, Geisser-Greenhouse, Huynh-Feldt, and the uncorrected, all use the same test statistic. For data analysis, UNIREP tests differ only by their respective degrees of freedom due to different degrees of freedom multipliers, which measure sphericity in the error covariance for the hypothesis variables. Muller and Stewart [6] provided detailed discussion of the basic theory for both the null and non-null cases. Earlier work detailed basic UNIREP theory. Box [7,8], Geisser and Greenhouse [9,10] and Huynh and Feldt [11] gave null results. Davies [12] and Muller and Barton [13,14] treated the non-null case.

Browne [15] evaluated the impact of using a pilot study to estimate the variance for a t-test. More generally, Taylor and Muller [16] demonstrated how to construct exact power confidence intervals for the general linear univariate model for a data-estimated variance and fixed means. The same authors also generalized the result to provide an exact confidence region around a power curve. Parallel results for the UNIREP setting would be equally useful. We generalize the methods in Taylor and Muller [16] to UNIREP tests for repeated measures. We use analytic and simulation results to demonstrate that the techniques allow computing approximate confidence intervals and regions for power with good accuracy for the UNIREP tests, based on an estimated covariance matrix and fixed means.

Existing results

A vector , (n×1), is lower case bold. A matrix, , is upper case bold with transpose ′, inverse −1 and generalized inverse −. Also, is an (n×1) vector of 1’s and is an (n×n) identity matrix. A diagonal matrix with (i,i) element z is written Dg (). The expected value, variance, and trace are E(), , and tr(), respectively. Throughout, Z∼χ2(ν,ω) indicates that Z has a noncentral chi-square distribution with ν degrees of freedom and noncentrality ω, while Z∼χ2(ν) indicates a central distribution. Similarly, Z∼F(ν1,ν2,ω) indicates X has a noncentral F distribution with ν1 numerator and ν2 denominator degrees of freedom, and noncentrality ω with cumulative distribution function F(ν1,ν2,ω). A central F is written Z∼F(ν1,ν2) with quantile q indicated . Writing indicates (p×1) is Gaussian with mean and covariance Σ (p×p). If (N×p) has independent rows and , then indicates follows a Wishart distribution with N degrees of freedom, covariance Σ, and noncentrality =E(′)E()Σ−1.

The general linear multivariate model,

assumes N independent rows and . In the model, is the fixed, known design matrix with 1≤rank(X)≤q, and contains fixed, unknown regression coefficients. For repeated measures ANOVA, one-group designs have rank ()=1, and two-group comparisons have rank ()=2. The associated general linear hypothesis is

suchthat definesthebetween-subject effects(rank a) while defines the within-subject effects (rank b). Requiring estimable Θ and full rank {,} ensure a testable hypothesis. Appropriate selections of the contrast matrices ( and ) and null matrix (Θ0) allows testing important one-degree-of-freedom parameters, such as the difference between two means, or a comparison of two trends.

For = (′)−′, unscaled noncentrality is = (Θ−Θ0)′−1(Θ−Θ0), scaled noncentrality is . Here Σ∗=′=Dg ()′ is the covariance matrix among the hypothesis variables, with ′=′=, and ={λ} the eigenvalues. Estimates are and , with ν=N−rank(), the error degrees of freedom. Furthermore, , and , with . The sum of squares hypothesis matrix is and the sum of squares error matrix is , which are independent of one another. The notation follows that in Muller and Stewart [6]. Additional notation is in Appendix A: Additional notation.

The univariate approach to repeated measures can be expressed in terms of the general linear multivariate model. The Box conservative (Box), the Geisser-Greenhouse (GG), the Huynh-Feldt (HF), and the uncorrected (Un) UNIREP tests use the same test statistic,

and a central F distribution to approximate the null distribution of T,

The sphericity parameter, ∊=tr2(Σ∗)/[btr(Σ∗2)], quantifies the spread of population eigenvalues and is used to discount the degrees of freedom. The term sphericity reflects the fact that uncorrelated Gaussian variables with equal variances in three dimensions have a spherical scattergram. The eigenvalues of Σ∗ are the variances of the (uncorrelated) principal components of the hypothesis response variables. Perfect sphericity requires ∊ = 1, which occurs with all eigenvalues equal. Minimal sphericity has ∊=1/b, which occurs with one nonzero eigenvalue. Other patterns of Σ∗ have 1/b<∊<1.

The Box conservative test uses the fixed, lower bound of ∊, while the uncorrected test uses the fixed, upper bound of ∊. With sphericity (∊=1), the uncorrected test is exact and uniformly most powerful (among similarly invariant tests). The Geisser-Greenhouse and Huynh-Feldt tests use the observed data to estimate ∊. The Geisser-Greenhouse estimator, , is the maximum likelihood (ML) estimator. The Huynh-Feldt estimator, was proposed as the ratio of two unbiased estimators. Their claim holds only for the special case of rank ()=1. Lecoutre [17] provided a more general form. In turn, Gribbin [18] and Chi et al. [19] described a rank-adjusted approximately unbiased estimator, , which applies to any general linear multivariate model. The rank-adjusted power approximation was shown through simulations to approximate observed mean power values as well as, or better than, the Huynh-Feldt power approximation (Chi et al. [19]). Only the rank-adjusted Huynh-Feldt estimator will be considered in the remainder of the paper.

Although the four UNIREP tests all use the same test statistic, they each use a different measure of sphericity, here indicated e. For data analysis, all four tests use a critical value . Here ν1=ab and ν2=bν. The Box test uses e=1/b, the GG test uses , HF uses , and the uncorrected test uses e=1. The p-value is then computed, for observed test statistic t, as p=1−F(t,ν1e,ν2e). In all cases . In turn, the p-values always have the reverse order, with the Box p-value being largest, and the uncorrected being smallest.

Muller et al. [5] showed that the distribution function of the UNIREP test statistic can be expressed exactly in terms of the distribution function of the sum of b positively and b negatively weighted independent chi-squares, namely y∼χ2(a,ω) and y∼χ2(ν),

Muller et al. [5] also reported accurate F approximations of the form

Here, y∗1∼χ2(ν∗1,ω∗), y∗2∼χ2(ν∗2), and . Parameters λ∗1, ν∗1, ω∗, λ∗2, and ν∗2 are defined in Appendix A: Additional notation. Power analysis involves {λ} and {ω}, with

the diagonal elements of the scaled noncentrality, ∗=′Dg ()−1=∗Dg ()−1.

Methods

Estimating approximate UNIREP power with estimated covariance and fixed means

By extending results in Muller et al. [5], the following lemma helps simplify the F approximations. Appendix B: Supporting lemmas and proofs contains all proofs.

Lemma 1

The constant in the critical value of the UNIREP test statistic approximation introduced by Muller et al. [5] is equal to 1,

Thus,

For known covariance and means, the power approximations for the Box, Geisser-Greenhouse, rank-adjusted Huynh-Feldt, and uncorrected tests are all of the form

Here, is equal to tr(Σ∗)/b with b equal to the rank of Σ∗, and . Table 1 contains values for e1 through e5 for the four UNIREP tests when , and . The expressions for ∊ and ∊ were derived using the properties described in Lemma B.1.

Table 1

Sphericity multipliers for UNIREP power approximations for fixed means

		Multipliers
Covariance	Test	e 1	e 2	e 3	e 4	e 5
Σ∗(known)	Un	1	1	∊ n	∊ d	∊ n
	HF	E∊~	E∊~	∊ n	∊ d	∊ n
	GG	E∊^	E∊^	∊ n	∊ d	∊ n
	Box	1/b	1/b	∊ n	∊ d	∊ n
Σ^_∗(estimated)	Un	1	1	∊~_n	∊^_d	∊~_n
	HF	∊~_r	∊~_r	∊~_n	∊^_d	∊~_n
	GG	∊^_d	∊^_d	∊~_n	∊^_d	∊~_n
	Box	1/b	1/b	∊~_n	∊^_d	∊~_n

In practice, some elements of may be estimated and hence random. The random elements imply random power values, as with estimated covariance and fixed means, , for , the unbiased restricted maximum likelihood (REML) estimator. A distinction must be carefully maintained between the estimation study and target study. The estimation study provides the covariance estimate and has sample size Nest, design matrix rank of rank (est), and νest=Nest−rank(est) degrees of freedom. The target study for which power is desired has sample size N, rank () and ν=N−rank() degrees of freedom.

The ML estimator from the Geisser-Greenhouse test, , is an obvious estimator for the target study’s ∊. For power analysis, a parallel estimator is available for ∊:

A better choice, given in the following lemma, uses a ratio of unbiased estimators. The result generalizes the rank-adjusted Huynh-Feldt estimator for data analysis. Appendix B: Supporting lemmas and proofs has derivations of moments as well as all proofs.

Lemma 2

For the non-null case, a ratio estimating ∊ in terms of correlated, but unbiased, estimators is [b]

The corresponding estimator for the null case is , the rank-adjusted Huynh-Feldt sphericity estimator.

For estimated covariance and fixed means, approximate estimated UNIREP power is

with , and e1 through e5 estimated if unknown (Table 1). Nearly every combination of , , , , 1 and 1/b was examined for each UNIREP test for the wide range of simulations discussed in Muller et al. [5]. The values chosen provided the most accurate results. In retrospect, they are natural choices as well.

Approximate UNIREP power confidence intervals

The solution to the UNIREP problem parallels the solution to the univariate problem in Taylor and Muller [16]. The methods apply to any general linear hypothesis, including one degree-of-freedom contrasts, such as pair-wise group comparisons and differences in linear trend between two groups. Tests giving scalar secondary parameters are also common for one-group designs and two-group comparisons. For known covariance and means, e5 is defined to be ∊ (Table 1), and the noncentrality in equation 9 is with and . For . Therefore, it follows that

For estimated covariance and fixed means, a ratio involving one biased and two unbiased estimators (Lemma B.2) for estimating λ∗1 may be written as

In parallel to the univariate setting, the distribution of can be approximated with a Satterthwaite approximation: with . Lower and upper tail probabilities, α and α, respectively, define the confidence coefficient, p=1−α−α. Also, and . Approximate confidence limits for the noncentrality may be calculated using the following:

Approximate lower and upper bounds are therefore and . The strict monotone dependence of the noncentral F function on the noncentrality ensures an approximate confidence interval for power. Lower and upper bounds on power are, with e1 through e4 defined in Table 1 for ,

and

Taylor and Muller [16] recommended one-sided power confidence intervals by noting that “the change from a one-sided to a two-sided confidence interval has little effect on the upper bound, but a large effect on the lower bound”. Muller and Fetterman [20] provided examples of a one-sided power confidence interval in the univariate case.

Approximate UNRIEP power confidence regions for power curves

The new methods allow calculating a confidence interval for a single power value. The logic of a proof in Taylor and Muller ([16] equations 2.1-2.13 and surrounding text) guarantees that accurate confidence regions are provided by the point-wise calculations. The proof may be sketched for the present setting as follows. Equations 14-21 establish the validity of the approximate confidence interval for a particular alternative hypothesis, as specified by the scalar constant tr(). The randomness in the noncentrality arises from a scalar random variable, , analogous to a variance. Equation 19 describes a single event with a specified probability. The inequality defining the event, and the associated probability, do not change for different values of the scalar constant tr(). The smooth and strictly monotone dependence of power on the noncentrality ensures the validity of equations 20-21. The proof is completed by noting that the monotonicity extends the simultaneity property to the power confidence region.

Figure 1 gives an example plot of approximate power confidence regions surrounding the predicted power curve for the rank-adjusted Huynh-Feldt test for ∊=0.720. Graphical representations such as Figure 1 help researchers accurately recognize the amount of uncertainty in their power calculation, and lead to better decisions about design.

Figure 1

Approximate 95% confidence region for predicted power of the rank-adjusted Huynh-Feldt test of interaction over tr() with and population for conditions described in Section ‘Simulation 1 with rank (one-group repeated measures ANOVA)’.

In some cases scientists prefer to consider sample size as a function of the pattern of mean differences. The theory already presented allows plotting sample size as a function of mean difference, albeit with a shift in algorithm. The power function must be numerically inverted to solve for the sample size desired. Taylor and Muller [16] outlined the steps of algorithm needed for the univariate case. Details are not presented here for the sake of brevity.

Results

Simulation overview

The accuracy of the new approximate confidence intervals isevaluatedfor a widerange of conditions. Appendix C: Simulation details contains more details of the simulations and examples. All simulations were conducted in SAS/IML (SAS 9.1, SAS Institute, 2003) using a version of LINMOD 3.4 modified to include the rank-adjusted Huynh-Feldt estimator and test. Predicted power values and approximate power confidence intervals were computed using a similarly modified version of POWERLIB 2.03. The modified versions of LINMOD and POWERLIB are available at .

Simulation 1 with rank (one-group repeated measures ANOVA)

The accuracy of the new approximate confidence intervals were evaluated for a completely within-subject design with p=9 repeated measures, N∈{10,20,40}, and q=rank()=1. Values for , contrast matrices , , and Θ0 were chosen to test a within-subject interaction for α=0.05. The model was chosen to ensure predicted power values for the Geisser-Greenhouse test of 0.20, 0.50, and 0.80, using the power approximation in Muller et al. [5]. Population covariance matrices were chosen to provide ∊∈{0.282,0.505,0.720,1.00}. The sphericity values were selected to cover a range of eigenvalue patterns (i.e., patterns of the principal component variances) arising from the structure of Σ∗. For example, if b=3, then gives ∊≈0.50. In turn, gives ∊=1/3≈0.33. Pseudo-random realizations of the error matrix, , were generated and tests were calculated. The observed mean power values for the four UNIREP tests were calculated and tabulated for 500,000 replications per condition.

For the conditions described above, additional pseudo-random realizations of the error matrix were generated using an estimating study with sample size, Nest, of 10 and rank of , rank (est), of 1 with 500,000 replications per condition for all four UNIREP tests. Corresponding estimated covariance matrices were calculated, as well as lower and upper bounds for power. Both one- and two-sided confidence intervals were evaluated with target coverages of 90% and 95%. The number of replications gave a standard error of estimated coverage probability less than or equal to 0.0003 for 1−α=0.95, and 0.0004 for 1−α=0.90, nearly guaranteeing 3 digits of accuracy. Only coverage of observed mean power values, and not predicted, was tabulated. The accuracy of the predicted power values, with respect to the observed, made it essentially redundant to consider both.

Only the worst case results for two-sided 95% confidence intervals are presented here. The worst cases occurred with the smallest sample size for the target study. Table 2 contains results for the Box conservative test with a target sample size of 10. For a wide range of sphericity values and target power values, the target 95% estimated coverage is consistently reached. The two cases in which the target coverage is not reached occur with large population sphericity and low power. Under these conditions, the Box conservative test would not be used in practice.

Table 2

Target 95% CI (two-sided) estimated coverage of simulated population power for

Test	∊	Power	Lower tail	Coverage	Upper tail
	0.282	0.123	1.1	97.8	1.1
		0.535	1.9	97.0	1.1
		0.930	1.7	97.3	1.0
	0.505	0.054	0.1	97.3	2.6
		0.266	0.5	97.0	2.5
Box		0.690	1.1	97.0	1.9
	0.720	0.052	0.4	94.1	5.5
		0.227	0.6	96.8	2.6
		0.569	1.4	97.0	1.6
	1	0.023	0.6	85.1	14.3
		0.117	0.5	96.0	3.5
		0.350	0.8	97.8	1.4
	0.282	0.155	3.1	94.7	2.2
		0.585	2.6	95.6	1.8
		0.942	1.8	96.6	1.6
	0.505	0.162	5.4	87.7	6.9
		0.520	3.8	90.6	5.6
GG		0.870	2.6	92.4	5.0
	0.720	0.203	2.4	92.3	5.3
		0.539	2.6	94.1	3.3
		0.856	3.3	94.2	2.5
	1	0.161	0.7	95.6	3.7
		0.438	1.4	97.0	1.6
		0.751	2.7	96.2	1.1
	0.282	0.166	3.8	93.5	2.7
		0.602	2.8	95.2	2.0
		0.946	1.9	96.3	1.8
	0.505	0.210	8.2	82.9	8.9
		0.592	4.7	88.5	6.8
HF		0.902	2.9	90.9	6.2
	0.720	0.271	3.6	90.9	5.5
		0.631	3.4	93.3	3.3
		0.904	4.0	93.6	2.4
	1	0.224	0.8	96.7	2.5
		0.531	1.8	97.1	1.1
		0.821	3.2	95.9	0.9

Standard error of coverage probability ×100≈0.0003×100.

Table 2 also contains coverage results for the Geisser-Greenhouse and the rank-adjusted Huynh-Feldt tests. The target 95% estimated coverage is consistently reached for extreme sphericity values for both tests. For midrange sphericity values, the coverage fell below the target coverage from 0.8% to 7.3% for the Geisser-Greenhouse, and 1.4% to 12.1% for the rank-adjusted Huynh-Feldt. Coverage accuracy improved as the estimated power increased. In practice, lower power values are of little concern. For target power of 0.80 for the Geisser-Greenhouse test, the largest deviation from the target 95% estimated coverage was 2.6% for the Geisser-Greenhouse test and 4.1% for the rank-adjusted Huynh-Feldt test. Both occurred for the population sphericity value of 0.505.

Only a spherical case is appropriate to consider for the uncorrected test because otherwise the test will have inflated test size. Simulation results in Table 3 show that the approximation for the uncorrected test (with sphericity) always reached the target estimated coverage for the uncorrected test. The conservative bias could be eliminated by using optimal maximum likelihood estimates for the common variance and covariance (Morrison [21]), rather than the unstructured covariance estimate. Additional small changes are needed, associated with degrees of freedom, and corresponding to making all choices of e1 through e5 equal to 1.

Table 3

Target 95% CI (Two-sided) estimated coverage of simulated population power for the uncorrected test

N	Population power	Lower tail	Coverage	Upper tail
10	0.238	0.5	97.5	2.0
	0.551	1.5	97.6	0.9
	0.835	3.2	96.1	0.7
20	0.215	0.8	97.3	1.9
	0.520	1.6	97.6	0.8
	0.814	3.1	96.2	0.7
40	0.207	0.9	97.2	1.9
	0.509	1.5	97.7	0.8
	0.806	3.0	96.3	0.7

Standard error of coverage probability ×100≈0.0003×100.

Although not presented here, in general, the accuracy of the coverage improved directly with increasing sample size, for all tests and conditions. The accuracy of the approximate confidence bounds for all four UNIREP tests also improved as the population sphericity increased.

Simulation 2 with rank ()

All of the simulations in the second example considered the condition of rank of greater than 1. The cases used p=5 repeated measures, N∈{16,32,48}, and q=rank()∈{2,4,8,16}, corresponding to a three-, five-, nine-, and seventeen-group comparison, respectively. Appropriate fixed matrices of regression coefficients, , contrast matrices, and , and Θ0 were chosen to test a within-subject interaction for a test size, α, of 0.05. The matrices were also chosen to ensure approximate target predicted power values for the rank-adjusted Huynh-Feldt test of 0.20, 0.50, and 0.80. Specific design matrices, , were defined. Population covariance matrices were chosen to provide specific population sphericity values, ∊∈{0.282,0.505,0.720,1.00}. Observed mean power values were simulated and tabulated in a similar manner to that described in section ‘Simulation 1 with rank (one-group repeated measures ANOVA)’.

Pseudo-random realizations of the error matrix were generated using an estimating study with sample size, Nest, of 16 and rank of , rank (est), of 4 with 500,000 replications per condition for all four UNIREP tests. Corresponding estimated covariance matrices were calculated, as well as lower and upper bounds for power using the methods presented in section ‘Approximate UNIREP power confidence intervals UNIREP power confidence intervals’. Approximate confidence interval coverage was defined as the proportion of the 500,000 simulated bound realizations that successfully covered the observed mean power values for each condition described above. Only coverage of observed mean power values, and not predicted, were tabulated. The accuracy of the predicted power values, with respect to the observed, made it essentially redundant to consider both. Both one- and two-sided confidence intervals were evaluated with target coverages of 90% and 95%.

In practical biomedical research, low power values are of little concern. Rarely will one have a power targeted below 0.70. Therefore, only the results for target power values of 0.80 will be presented and discussed. Power confidence interval coverage converged to the target coverage as sample size increased. Only the worst case results for two-sided 95% confidence intervals are presented here. The worst cases occurred with the smallest sample size for the target study, for a variety of population sphericity values and estimated population powers.

In Table 4, the observed mean population powers are presented for the four UNIREP tests for the population sphericity values and ranks of considered for target rank-adjusted Huynh-Feldt power of 0.80 and sample size of 16 or 48. In general, as the population sphericity increased and rank of increased, the observed mean power values for the Box conservative, the Geisser-Greenhouse, and the rank-adjusted Huynh-Feldt tests decreased. Only the Box conservative had severely biased power values as the population sphericity increased.

Table 4

Simulated population power for target power and rank ( )

N				Simulated population power
		∊=0.282			∊=0.505
	q	Box	GG	HF	Box	GG	HF
16	2	0.779	0.811	0.817	0.561	0.778	0.809
	4	0.763	0.797	0.805	0.510	0.762	0.802
	8	0.753	0.787	0.799	0.455	0.736	0.796
		∊=0.720			∊=1.000
	q	Box	GG	HF	Box	GG	HF	UN
16	2	0.457	0.760	0.805	0.399	0.748	0.790	0.799
	4	0.355	0.740	0.801	0.255	0.724	0.787	0.801
	8	0.267	0.695	0.795	0.138	0.655	0.775	0.800
		∊=0.282			∊=0.505
	q	Box	GG	HF	Box	GG	HF
48	2	0.803	0.843	0.845	0.586	0.802	0.813
	4	0.773	0.812	0.815	0.552	0.794	0.805
	8	0.766	0.800	0.803	0.544	0.787	0.799
	16	0.762	0.796	0.799	0.522	0.780	0.795
		∊=0.720			∊=1.000
	q	Box	GG	HF	Box	GG	HF	UN
48	2	0.500	0.792	0.806	0.455	0.785	0.797	0.800
	4	0.427	0.787	0.802	0.346	0.781	0.796	0.800
	8	0.402	0.781	0.799	0.280	0.778	0.797	0.801
	16	0.359	0.775	0.799	0.221	0.770	0.795	0.801

Standard error of estimated power ≈0.0006.

In Table 5, the proportion of simulations in which the estimated confidence interval successfully covered the observed mean population power values for each test is shown. The results are based on using an estimating study with sample size, Nest, of 16 and rank of , rank (est), of 4. In general, the approximate power confidence intervals nearly always reached the target 95% coverage for the Box conservative test. The coverage became more conservative as rank of decreased. Similarly, the coverage became more conservative for the Geisser-Greenhouse and rank-adjusted Huynh-Feldt tests as rank of decreased. The Geisser-Greenhouse and rank-adjusted Huynh-Feldt tests performed adequately in all cases except for the midrange population sphericity value, ∊=0.505. The largest deviation from the target 95% estimated coverage was 13.6% and 16.0% for the Geisser-Greenhouse and rank-adjusted Huynh-Feldt tests, respectively, which occurred for ∊=0.505 and rank of equal to 8. The approximate power confidence intervals for the uncorrected test reached the target coverage value for every case considered in which the uncorrected test would be used.

Table 5

Target 95% CI (two-sided) estimated coverage of simulated population power for target power rank N and rank ( )

	∊ = 0.282					∊ = 0.505
q	Box	GG	HF		Box	GG	HF
2	97.8	97.2	97.0		97.5	93.4	92.3
4	93.7	92.0	91.6		95.6	86.8	85.0
8	90.9	87.9	87.2		94.8	81.4	79.0
	∊=0.720					∊=1.000
	Box	GG	HF		Box	GG	HF	UN
2	97.6	95.4	94.9		97.4	95.3	95.5	95.8
4	97.5	93.6	92.9		97.6	96.8	97.0	97.4
8	96.9	90.6	89.8		97.0	96.1	96.9	97.4

Standard error of coverage probability × 100 ≈ 0.0003 × 100.

Although not presented here, in general, as sample size increased the conservative coverage values observed for the Box conservative and the uncorrected tests slowly converged to the target coverage value. This trend was observed for the conservative coverage values with the extreme population sphericity values for the Geisser-Greenhouse and the rank-adjusted Huynh-Feldt tests as well. The same is true of the liberal coverage values observed for the midrange population sphericity values for the Geisser-Greenhouse and the rank-adjusted Huynh-Feldt tests. Similar results were obtained for the target 90% two-sided confidence interval coverage as well as the 95% and 90% one-sided confidence intervals coverage.

The estimated coverages of these tabulated observed mean power values for each test were simulated for population sphericity values of 0.282, 0.505, 0.720, and 1.00. In Table 6, the worst case results from these simulations, which occurred for population sphericity 0.505, are presented. Approximate confidence intervals were simulated for 500,000 replications per condition (standard error of estimated coverage probability less than or equal to 0.0003 for 1−α=0.95, and 0.0004 for 1−α=0.90). The estimating studies use sample sizes, Nest, of 16, 32, and 48, and ranks of est of 2, 4, and 8.

Table 6

Target 95% CI (two-sided) estimated coverage of simulated population power for , target power and rank

		Box coverage			GG coverage			HF coverage
		rank ( X est )			rank ( X est )			rank ( X est )
N est	q	2	4	8	2	4	8	2	4	8
16	2	97.5	97.2	97.4	94.1	94.3	94.9	92.2	92.1	93.6
	4	94.8	94.8	95.3	87.3	87.5	88.9	85.9	86.2	87.4
	8	92.6	92.7	93.4	83.2	83.4	86.0	81.4	82.0	84.3
	16	92.0	92.3	93.7	82.3	82.5	85.9	80.5	80.7	83.2
32	2	97.3	97.3	97.3	93.3	93.3	93.4	92.4	92.5	92.6
	4	93.8	94.1	94.3	85.7	85.2	85.8	85.0	84.8	84.5
	8	91.6	91.8	91.4	81.3	81.5	82.4	79.5	80.0	80.6
	16	91.5	91.7	91.7	79.4	78.9	80.0	79.2	79.1	79.5
64	2	97.2	97.2	97.4	93.6	93.7	93.5	92.6	92.5	92.8
	4	94.4	94.6	94.8	84.5	85.0	84.7	84.4	85.0	85.2
	8	91.7	91.5	91.8	79.6	80.1	80.4	78.9	79.2	79.6
	16	90.9	90.9	91.0	78.5	78.4	78.7	78.9	78.4	78.7

Estimation Study: (16, 32, 48) and rank () (2,4,8). Standard error of coverage probability × 100 ≈ 0.0003 × 100.

In general, for population sphericity values of 0.282 and 0.505, the approximate power confidence interval coverage for the Box conservative test converged to the target coverage value as rank of est increased, and thus νest decreased. Coverage decreased as rank of from the target study increased. For larger rank of , the approximate power confidence interval coverage fell short of the target coverage in several instances. No clear trend was apparent as Nest increased. The Box conservative test would not be used for larger population sphericity values. However, the realization that the target coverage was reached in nearly every case considered for the larger population sphericity values is worth mentioning.

The approximate power confidence interval coverages for both the Geisser-Greenhouse and rank-adjusted Huynh-Feldt tests seem to have converged to the target coverage value as rank of est increased, and thus νest decreased, except in cases of sphericity. Such cases have little practical importance since exact results are available if sphericity is valid. Coverage decreased as rank of from the target study increased. As observed in previous simulations, the approximate power confidence interval coverages for both the Geisser-Greenhouse and rank-adjusted Huynh-Feldt tests fell short of the target coverage to varying degrees in nearly every case considered for midrange population sphericity values. This outcome was also observed for larger rank of from the target study for population sphericity of 0.282. The approximate power confidence interval coverage for the uncorrected tests reached the target coverage value in every case except for large νest and small rank of from the target study. The approximate power confidence interval coverage increased as the ranks of for both the target and estimating studies increased and as Nest decreased.

The slow convergence of the approximate power confidence interval coverage to the target coverage may be due, in part, to use of and in the approximate power confidence interval equation. These estimators of the sphericity parameter are ratios of unbiased estimators for the non-null and null cases, respectively. The variances of these estimators are much larger than the variances for and . The larger variances may account for the slow convergence to the population power as the target and estimating study sample sizes and degrees of freedom increase. Further simulations may be needed to confirm this reasoning.

Alternate approximations

In attempts to develop even better confidence bound estimates, additional approximations were developed and evaluated. One approach approximated the distribution of with an F. Using the methods presented in Kim et al. [22], the numerator of was approximated with a weighted noncentral chi-square, while the denominator was approximated with a weighted central chi-square. Two concerns arose. First, the denominator is not necessarily a central quadratic. The 2tr(/a) component makes the denominator more of a shifted central quadratic. Second, the Kim et al. [22] result requires that the components of the numerator and denominator be mutually independent, which does not hold. Simulations demonstrated that the approximation was inaccurate in small samples.

Alternative approximations were developed and evaluated. The alternatives matched only the numerator to a weighted noncentral chi-square or to a weighted central chi-square with the denominator a constant equal to E. All were less accurate than the approximation presented here.

Discussion

Even for small sample sizes, the proposed power confidence intervals attain very accurate coverage probabilities for the Box conservative test in all cases and for the uncorrected test with ∊=1 (the only case for which it should be used). The result is also true for the extreme population sphericity values for the Geisser-Greenhouse and rank-adjusted Huynh-Feldt tests. For midrange population sphericity values, the coverage probabilities of the approximate power confidence intervals for the latter two tests often fell somewhat short of the various target coverage values considered. Coverage probabilities improve as sample size increases. Accuracy is better for higher target power values than for lower, which makes the results useful in practice. One-sided confidence intervals are recommended for lower bounds on power.

The techniques also allow plotting power confidence regions around an estimated power curve (Figure 1). The resulting plots have been well received by researchers.

Conclusions

Good statistical practice requires associating a credible measure of uncertainty with any parameter estimate. We described and evaluated new methods to meet the need for UNIREP power estimates based on an estimated covariance with fixed means. Across a large range of conditions, the methods provide reasonably accurate coverage for all four UNIREP tests.

Appendix

Appendix A: Additional notation

The additional notation comes from Muller et al. [5], who showed that if , , and , then

They used the parameters (assumed known) to approximate the UNIREP test statistic with a noncentral F distribution, as presented in equation 6.

Appendix B: Supporting lemmas and proofs

Lemma B.1

A.1-A.5 imply .

Proof of lemma B.1

Lemma B.2

The first moments of , , , and are known.

Proof of lemma B.2

Following Wishart [23], , such that ν=N−rank(). For Wishart while here 〈Σ〉=σ and ρ=σ/(σσ)1/2. Wishart [23] said, with emphasis not in the original, “ …moment coefficients are in all cases except the first calculated about the mean of the sample …”. Here μ(n) indicates the expression in equation n at the end of Wishart [23], . Thus

With , . In turn . Also:

Hence

Here , with . Thus,

Finally, has E ()=νΣ∗ (Muller and Stewart [6], Theorem 10.10). Hence E [tr()]=tr[E()]=tr[E()]=tr[(νΣ∗)]=νtr(∗) and

Proof of lemma 1

Substituting equivalent terms from equations A.1-A.5 into (λ∗2abν∗2)/(λ∗1ν∗1bν) allows combining like terms and simplifying to give

If , then and with y∗1∼χ2(ν∗1,ω∗) and y∗2∼χ2(ν∗2) as described in Muller et al. [5]. In turn,

Pr{(y∗1/ν∗1)/(y∗2/ν∗2)≤tλ∗2abν∗2/(λ∗1ν∗1bν)}=F(t;ν∗1,ν∗2,ω∗).

Proof of lemma 2

Using Lemma B.2, unbiased estimators for tr2(Σ∗) and are and , as introduced in Gribbin [18] and Chi et al. [19]. As shown in Lemma B.2, and are unbiased estimators for and , respectively. Thus,

In the null case = and reduces to the rank-adjusted sphericity estimator, .

Appendix C: Simulation details

Covariance conditions

Covariance conditions 5-8 from Table III of Coffey and Muller [24] were used for each example described below: Σ∗=Dg() for j∈{1,2,3,4}, with

Thus, ∊∈{0.28,0.51,0.72,1.00}. Given Σ∗=Dg(), it follows that Σ=∗′.

CLAHE mammography example

Computer scientists developed the Contrast Limited Adaptive Histogram Equalization (CLAHE) algorithm to improve contrast in medical images. Independent observers considered 3×3=9 Clip ×Region combinations. Region denotes the size of the image (pixels2) at which contrasts are controlled and Clip level limits the maximum contrast adjustment. In the multivariate model =, while within-person factors Clip and Region gave , (N×9). Also , (1×9), contained mean log10(contrast) for the unprocessed condition minus the mean for each of the nine combinations of Clip and Region (βcr=μunprocessed−μcr). If c contains orthonormal linear and quadratic trends for log2(Clip) ∈ {1,2,4}, and r does the same for log2(Region) ∈ {1,3,5}, then the 9 × 4 within-persons contrast matrix, cr is

With L the linear and Q the quadratic trends for interaction components being tested, .

All four covariance patterns were factorially combined with N∈{10,20,40}. The multivariate test considered with α=0.05, and β the scaling factor for corresponding to approximate target power P∈{0.20,0.50,0.80} for the Geisser-Greenhouse approximation using methods in Muller et al. [5]. The conditions in the example were used in section ‘Simulation 1 with rank (one-group repeated measures ANOVA)’. More details of the example are in Muller et al. [5].

Test of interaction with rank () Example

All cases used 5 repeated measures, N∈{16,32,48}, and rank ()∈{2,4,8,16}. For obvious reasons, a rank of equal to 16 was not considered for the smallest sample size. All four covariance patterns were factorially combined with the sample sizes and ranks . In the multivariate model, =⊗, such that repn = N/q, and ⊗ is a Kronecker product. If

then =β·, such that β was the scaling factor giving approximate target power P∈{0.20,0.50,0.80}, for the rank-adjusted Huynh-Feldt power approximation. The within-subject contrast, , (5×4), contained orthonormal linear, quadratic, cubic and quartic trends:

Each row of the between-subject contrast, , a (q−1×q) orthonormal matrix, contained one of the (q−1) trends. The contrasts define a test of interaction of between- and within-subject trends. Without loss of generality, we assumed Θ0=. A test size, α, of 0.05 was used.

Computational methods

All power computations were conducted in SAS/IML (SAS 9.1, SAS Institute, 2003). Free software LINMOD 3.4 was used for all data analysis and includes new methods. Free software POWERLIB 2.1 in Johnson et al. [25] was used for all power analysis and includes the new methods. Both are available at . UNIREP power is also available in GLIMMPSE, a free web-browser based program with a graphical user interface aimed at health scientists (). The next version of GLIMMPSE is expected to implement the confidence interval methods.

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

The majority of the work comes directly from the dissertation of MJG at UNC Chapel Hill, in partial fulfillment of the requirements of the DrPH, under the direction of KEM. Dr. KEM and Dr. Y-YC assisted in revising and compressing the notation, text and tables. Dr. PWS was a member of the dissertation committee and assisted in formulating the proof of simultaneous coverage of the confidence regions and in final editing. All authors read and approved the final manuscript.

Pre-publication history

The pre-publication history for this paper can be accessed here:

http://www.biomedcentral.com/1471-2288/13/57/prepub